AD-A213  672 


f\to  3  i  ny.io '6-s 


a 


Final  Report 

Experimental  and  Theoretical  Studies  of  Wave  Propagation 
in  Granular,  Rock  and  Porous  Media 


by 

Martin  H.  Sadd 
Arun  Shukla 


Prepared  for 

The  U.S.  Army  Research  Office 
Under  Contract  No.  DAAL03-86-K-0125 


SLECTE 
OCT  2  5  1939 


Department  of  Mechanical  Engineering  &  Applied  Mechanics 

University  of  Rhode  Island 
Kingston,  RI  02881 


August  1989 


-f 

Jl 


MASTER  COPY 


FOR  REPRODUCTION  PURPOSES 


UNCLASSIFIED 

C’jRirv  classificat; 


ia.  report  security  CLASSIFICATION 


REPORT  DOCUMENTATION  PAGE 

lib.  RESTRICTIVE  MARKINGS 


[  2*.  SECURITY  CLASSIFICATION  AUTHORITY 


|  2b.  DECLASSIFICATION  /  DOWNGRADING  SCHEDULE 


I  4.  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 


6a.  NAME  OF  PERFORMING  ORGANIZATION 

University  of  Rhode  Island 


6b.  OFFICE  SYMBOL 
(If  apfilicaoie) 


6c  AOORESS  (Cty,  Start,  and  21?  Cod*) 

Wales  Hall 
Kingston,  RI  02881 

8a.  NAME  OF  FUNDING /SPONSORING  I  3b.  OFFICE  SYMBOL 

ORGANIZATION  Qf  apfili cable) 

U.  S.  Army  Research  Office 

8c  AOORESS  (Cty,  Start,  and  ZIP  Code) 

P.  0.  Box  12211 

Research  Triangle  Park,  NC  27709-2211 


3.  DISTRIBUTION  /  AVAILABILITY  OF  REPORT 

Approved  for  public  release; 
distribution  unlimited. 

5.  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 

S.22X2./Q  _ 

7a.  NAME  OF  MONITORING  ORGANIZATION 

U.  S.  Army  Research  Office 

7b.  AOORESS  (Cry,  Start,  and  ZIP  Code) 

P.  0.  Box  12211 

Research  Triangle  Park,  NC  27709-2211 

9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 

J)  k-0  _ 

10.  SOURCE  OF  FUNDING  NUMBERS 


PROGRAM 
ELEMENT  NO. 


PROJECT 

NO. 


WORK  UNIT 
ACCESSION  NO 


11.  TITLE  (Include  Security  Classification) 

Experimental  and  Theoretical  Studies  of  Wave  Propagation  in  Granular,  Rock  and  Porous 
Media. 


13b.  TIME  COVERED 
FROM  7/1/86  TO  6/ 


14.  DATE  OF  REPORT  (Year,  Month,  Day)  15.  PAGE  COUNT 

19  1989,  August,  31 _ 132 


12.  PERSONAL  AUTHOR(S) 

Martin  H.  Sadd,  and  Arun  Shukla 

13a.  TYPE  OF  REPORT  13b.  TIME  COVERED  14.  DATE  OF  REPORT  (far.  Month,  Day)  15.  PAGE  COUNT 

Final  Report _  from  7/1/86  TO  6/30/89I  1989,  August,  31  132 _ 

16.  SUPPLEMENTARY  NOTATION 

.  The  view,  opinions  and/or  findings  contained  in  this  report  are  those 

of  the  author (s) .and  should  not  be  construed  as,  an  official  Department  of  the  Army  position 

nnl  1  f'L  -QT  ilpn-iinn.  nnTp^q  an  Hpt-j  gnnrorl  Hy  nfhgr-  flflCUafiaLa&iaa 

17 _ _ COSATI  COOES  18.  SU8JECT  TERMS  (Continue  on  Havana  if  necessary  and  identify  by  block  number) 

FIELD _ GROUP _ SUB-GROUP  _  ^»Wave  Propagation,  Granular  Media,  Microstructural  Media, 

_ S  Dynamic  Photomechanics,  Distinct  Element  Method,  Rock 

Mechanics  &  Porous  Media.  _ 

1  ^ABSTRACT  (Continue  on  reverse  if  necessary  and  identify  by  block  number)  S 

An  experimental  and  theoretical  study  of  wave .propagation  in  granular’,  rock  and 
porous  media  is  described.  The  wave  propagation  is  a“~result_o_f  explosive  loadings 
.of  50-100  duration.  Specific  wave  propagational  variables  include  the  wave 
speed,  intrer-granular  contact  forces  (local  wave  amplitude)  and  wave  spreading 
geometrie/s.  The  relationship  of  the  media  microstructure  to  these  propagational 
variables  has  been  studied  for  a  variety  of  granular  and  porous  materials.  \ 


20  C'  -RIBUTION/ AVAILABILITY  OF  ABSTRACT  21.  ABSTRACT  SECURITY  CLASSIFICATION 

□  UNCLASSIFIEDAJNUMITED  □  SAME  AS  RPT.  □  OTIC  USERS _ Unclassif  ied _ 

lla.  NAME  OF  RESPONSIBLE  INDIVIDUAL  22b.  TELEPHONE  (Include  Area  Code)  I  22c.  OFFICE  SYMBOL 


OO  FORM  1473,  84  MAR 


83  APR  edition  may  b«  used  until  exhausted. 
All  other  editions  are  obsolete. 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 

UNCLASSIFIED 


Table  of  Contents 


Report  Documentation  Page 
Table  of  Contents 

1.  Statement  of  Problem 

2.  Summary  of  Results 

3.  List  of  Publications 

4.  List  of  Participating  Scientific  Personnel 

5.  Appendices:  Reprints  of  Major  Publications 


1.  Statement  of  Problem 


The  completed  research  to  be  described  in  this  report  is  concerned  with  the 
propagation  of  waves  in  geological  materials.  Since  such  materials  are  made  up  from 
sands,  clays  and  rocks,  they  contain  granular  and  porous  microstructures.  These  local 
microstructures  or  fabric  (as  it  is  sometimes  referred  to)  are  a  function  of  the  geometry 
of  the  material  distribution,  and  produce  non-continuous  fields  of  mass  density,  stress, 
strain,  displacement,  etc.  These  types  of  materials  are  therefore  difficult  to  model  using 
classical  continuum  mechanics  such  as  elasticity,  plasticity  or  viscoelasticity  theories. 
One  of  the  key  features  of  this  research  was  to  investigate  the  effect  of  this 
microstructure  on  the  propagation  of  waves  produced  by  explosive  loadings  of  short 
duration. 

In  many  cases,  geological  materials  will  be  composed  of  sands,  gravel  and/or 
crushed  rock,  and  can  therefore  be  categorized  as  a  granulai  num.  Such  a  medium 
can  be  characterized,  for  modeling  purposes,  as  a  collection  of  distinct  particles  which 
can  displace  independently  from  one  another  and  interact  only  through  contact 
mechanisms.  This  type  of  media  transmits  mechanical  loadings  through  discrete  paths 
as  determined  by  the  geometry  of  the  granular  packing.  For  the  case  of  porous  media 
containing  various  distributions  of  open  pores,  a  similar  direction-oriented  behavioral 
response  will  also  occur. 

In  regard  to  wave  propagation,  our  interest  was  involved  with  the  transient 
dynamic  behavior  of  such  a  structured  medium  when  it  is  subjected  to  explosive 
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loadings  of  50-100|is  duration.  Specific  variables  of  interest  included: 


•  wave  speed 

•  inter-granular  contact  forces  (local  wave  amplitude) 

•  wave  spreading  geometry 

The  goal  of  the  research  was  to  relate  how  the  specific  microstructure  of  the  medium 
effects  these  wave  propagational  variables. 

In  order  to  conduct  the  investigation,  both  theoretical  and  experimental  methods 
of  study  were  employed.  The  theoretical  work  focused  on  two  general  methods  which 
included  the  use  of  a  new  microstructural  continuum  theory  for  modeling,  and  a 
computational  study  using  the  so-called  distinct  element  method  to  model  the  dynamic 
behavior  of  granular  assemblies.  The  microstructural  continuum  theory  which  was  used 
was  the  distributed  body  theory  originally  developed  by  Goodman  and  Cowin  (see 
reference  29  in  paper  A  in  the  Appendix).  This  particular  theory  assumes  that  the 
medium  is  distributed  in  space  by  an  independent  kinematical  function  called  the 
volume  distribution  function,  and  thus  this  theory  allows  the  medium  to  contain  voids 
and  other  microstructures.  In  regard  to  the  numerical  modeling,  the  distinct  element 
method  was  used  to  determine  the  motion  of  each  granule  in  particular  model 
assemblies.  This  computational  scheme  assumes  that  each  granule  may  be  modeled  by 
rigid  body  dynamics  with  particle  interactions  having  stiffness  and  damping  properties. 
The  predicted  movements  of  the  various  granules  are  then  the  result  of  the  propagation 
through  the  medium  of  disturbances  which  originated  at  the  boundary  loading  points. 

The  experimental  efforts  employed  the  principles  of  dynamic  photoelasticity  along 
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with  high  speed  photography  to  get  the  complete  stress  field  information  due  to 
explosive  loading  in  various  granular  assemblies.  This  technique  gave  complete  demils 
about  the  wave  velocities,  dynamic  contact  stresses  and  wave  spreading  during  the 
particle  contact  interaction  period.  Assemblies  containing  particles  of  different  size  and 
arrangements  were  used  in  the  experiments.  A  detailed  investigation  was  conducted  to 
study  the  influence  of  local  microstructure  on  the  wave  propagation  phenomenon.  The 
experimental  study  also  made  use  of  electrical  resistance  strain  gages  to  study  wave 
propagation  in  real  earth  materials.  A  preliminary  investigation  was  conducted  to  study 
dynamic  load  transfer  in  granular  rock  media. 
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2.  Summary  of  Results 


This  section  briefly  describes  the  major  results  of  the  research  program. 

References  will  be  made  to  the  reprints  of  our  major  papers  which  appear  in  the 

Appendix.  Details  of  the  particular  results  can  be  found  in  each  of  the  papers. 

Wave  Propagation  in  Distributed  Bodies  Wave  propagation  studies  have  been 
conducted  based  upon  a  distributed  body  model  of  geological  media.  The 
distributed  body  model  employs  an  independent  kinematical  volume  distribution 
function  which  describes  the  way  the  material  is  distributed  in  space,  and  thus 
allows  the  theory  to  generate  porous  and  granular  microstructures.  The  theory 
uncouples  the  mass  density  of  the  granules  from  the  mass  density  of  the  entire 
material,  and  allows  compressibility  due  to  both  granule  compressibility  and  void 
compaction.  A  one-dimensional  theory  of  propagating  singular  surface 
acceleration  waves  has  been  developed  for  particular  volume  distribution 
functions  which  have  application  to  geological  materials.  Three  different  volume 
distribution  functions  were  developed  producing  periodic,  exponential  and  a 
periodic-exponential  material  microstructures.  A  general  computer  code  was 
developed  to  calculate  the  results  predicted  from  the  model  for  a  variety  of 
constitutive  and  micro  structural  model  parameters.  Results  for  the  wave  speed 
and  amplitude  behavior  as  a  function  of  these  model  parameters  are  given  in 
Figures  2-8  in  Paper  A  in  the  Appendix.  Specific  relationships  between  the 
microstructure  and  the  wave  speed  and  attenuation  have  been  determined,  and 
these  theoretical  results  are  in  general  agreement  with  measured  results. 


Distinct  Element  Modeling  of  Wave  Propagation  in  Granular  Materials  The 
distinct  element  model  is  a  numerical  scheme  which  uses  Newtonian  rigid-body 
mechanics  to  model  the  translational  and  rotational  motion  of  each  disk  in  a 
granular  assembly.  The  method  incorporates  prescribed  stiffness  and  damping 
at  the  contacts  between  each  of  the  granules,  and  these  model  parameters  were 
determined  from  calibration  tests  in  the  experimental  segment  of  the  research 
program.  Ultimately  the  distinct  element  method  establishes  an  explicit  time¬ 
stepping  scheme  that  enables  the  calculation  of  the  inter-granule  contact  force 
between  all  granules  at  each  of  the  selected  time  steps.  The  method  is 
numerically  efficient  so  that  large  numbers  of  granules  may  be  analyzed  in  a 
given  problem.  The  basic  method  is  outlined  in  Paper  B  in  the  Appendix. 
Comparisons  of  the  numerical  results  with  those  from  dynamic  photoelasticity  are 
shown  in  Figures  4  and  5  in  Paper  B.  These  comparisons  indicate  reasonably 
good  agreement  between  theoretical  and  experimental  results  for  the  assemblies 
which  were  studied.  Inter-granular  contact  force  distributions  indicated  the 
dependence  of  microstructure  on  the  wave  propagational  characteristics,  and  it 
was  discovered  that  the  microstructural  measure  of  the  branch  angle  between 
local  granules  is  an  effective  variable  to  use  in  establishing  the  connection 
between  microstructure  and  wave  propagation. 

The  Effect  of  Voids  and  Inclusions  on  Wave  Propagation  in  Granular 
Materials  Theoretical  and  experimental  studies  have  been  conducted  on  wave 
propagation  in  granular  materials  containing  local  discontinuities  of  voids  and 
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inclusions.  The  granular  medium  was  simulated  by  specific  assemblies  of 
circular  disks,  and  the  voids  were  created  by  removing  particular  disks  from  the 
assembly  while  inclusions  were  created  by  replacing  certain  disks  with  ones  of 
a  higher  impedance  material.  The  theoretical  modeling  employed  the  distinct 
element  method,  and  the  experimental  study  used  the  technique  of  dynamic 
photoelasticity.  Comparisons  were  made  between  the  computational  and 
experimental  data  on  the  inter-granular  contact  forces  around  each  void  or 
inclusion,  and  these  comparisons  are  shown  in  Figures  7-11  in  Paper  C  in  the 
Appendix.  Both  voids  and  inclusions  produce  local  wave  scattering  through 
various  reflection  mechanisms,  and  the  results  seem  to  indicate  that  the  inclusions 
produce  higher  local  wave  attenuation. 

Experimental-Numerical  Hybrid  Technique  of  Load  Transfer  Coefficients  for 
Wave  Propagation  Predictions  in  Granular  Media  Experimental  studies  on 
dynamic  load  transfer  in  granular  media  have  been  conducted  through  the  use 
of  dynamic  photoelasticity.  The  experimental  data  collected  allows  the 
determination  of  the  time  dependent  inter-granular  contact  loadings  between  the 
granules.  This  leads  to  the  calculation  of  load  transfer  coefficients,  i.e.  the  ratio 
of  the  maximum  output  contact  load  to  the  maximum  input  contact  load,  for 
various  packing  geometries.  These  coefficients  were  then  used  along  with  the 
principle  of  superposition  to  predict  the  peak  inter-granular  contact  loads  in 
several  model  granular  assemblies.  Results  of  the  numerical  hybrid  scheme  were 
then  compared  with  experimental  data  for  the  assemblies  investigated,  see  Figures 
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12-15  in  Paper  D  in  the  Appendix. 

Influence  of  Local  Microstructure  on  Wave  Propagation  Phenomenon  A 
detailed  experimental  study  was  conducted  to  evaluate  the  effect  of  different 
microstructures  or  fabric  of  granular  media  on  the  wave  propagation  phenomenon. 
The  granular  media  was  simulated  by  circular  disks  made  of  photoelastic 
materials.  Attention  was  focussed  on  the  load  transfer  paths,  wave  velocities, 
wave  attenuation  and  the  dynamic  stresses  which  are  generated  at  the  contacts 
due  to  the  passage  of  stress  waves.  The  details  of  the  results  are  included  in 
Paper  E  in  the  Appendix.  Figures  2  through  7  contain  photographs  showing  the 
full  field  description  of  the  wave  propagation  process  in  different  assemblies. 
The  paper  discusses  the  primary  and  secondary  load  transfer  paths  in  different 
assemblies  and  clearly  demonstrates  the  influence  of  the  microstructure  on  the 
load  transfer  process. 

Angular  Dependence  of  Dynamic  Load  fransfer  Process  An  experimental 
investigation  was  conducted  to  evaluate  the  dependance  of  packing  geometry  on 
the  dynamic  load  transfer  in  two  dimensional  granular  chains.  The  results  from 
these  experiments  are  shown  in  Paper  F  in  the  Appendix.  It  was  shown,  that 
in  two  dimensional  chains,  rapid  attenuation  of  load  transfer  occurs  as  the  branch 
angle  increases  from  0  to  90  degrees.  It  was  also  observed  that  the  wavelength 
of  the  loading  pulse  increases  with  the  branch  angle.  These  results  are  illustrated 
in  the  Figures  1 1  and  12  of  Paper  F. 
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Wave  Propagation  in  Porous  Media  As  a  Function  of  Fluid  Saturation 
Dynamic  photoelasticity  was  used  to  study  wave  propagation  in  a  porous  media 
as  a  function  of  fluid  saturation.  The  porous  media  was  modeled  as  a 
continuous  solid  containing  particular  arrays  of  holes  or  voids.  The  study 
investigated  the  wave  propagation  phenomenon  from  a  microscopic  point  of  view 
by  going  into  the  details  of  the  •  eometric  nature  of  the  porous  structure.  The 
details  of  this  study  are  given  in  Paper  G  in  the  Appendix.  The  results  show 
the  dependence  of  wave  velocity  and  attenuation  on  the  porosity  as  well  as  the 
microstructural  arrangement  of  the  pores. 

Wave  Propagation  and  Energy  Transfer  Across  Contacts  B etv.  ee n  Large  Bodies 
This  study  experimentally  investigated  the  formation  of  dynamic  contacts  between 
two  bodies  which  are  much  larger  than  the  wave  length  of  the  loading  pulse. 
Results  from  this  study  are  given  in  Paper  H  in  the  Appendix.  Dynamic 
photographs  provided  full  field  information  of  diffraction,  reflection  and 
transmission  process  as  a  function  of  time.  The  results  show  that  the  individual 
wave  types  interact  with  the  contact  region  in  specific  ways  as  determined  by 
the  reflection  and  refraction  laws.  It  must  be  mentioned  here  that  only 
qualitative  treatment  of  the  data  was  possible  in  this  study.  Initial  attempts  to 
quantitatively  evaluate  fringe  pattern  data  were  not  successful. 

Dynamic  Load  Transfer  in  Virgin  and  Damaged  Rock  Media  Electrical 
resistance  strain  gages  were  used  to  study  dynamic  load  transfer  in  a  single  chain 
assembly  of  disks  fabricated  from  four  different  grades  of  white  Vermont  Marble. 
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The  study  was  preliminary  in  nature  and  the  results  are  given  in  Paper  I  in  the 
Appendix.  The  results  show  the  dependence  of  wave  velocity  on  the 
microstructure  of  the  ^ock.  The  effect  of  prior  damage  in  rock  grains  on  wave 
propagation  phenomenon  was  also  investigated. 
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Sadd  and  Hossain 


Paper  A 


Wave  Propagation  in  Distributed  Bodies  ... 


WAVE  PROPAGATION  IN  DISTRIBUTED  BODIES  WITH 
APPLICATIONS  TO  DYNAMIC  SOIL  BEHAVIOR 

MARTIN  H.  SADD  and  MOHAMMAD  HOSSAIN 
Mechanical  Engineering  and  Applied  Mechanics 
University  of  Rhode  Island 

Kingston,  RI  U  2  jii  -  0  T 

ABSTRACT 

A  study  is  presented  which  models  wave  propagation  through  materials  with  microstructure.  The 
specific  microstructure  of  interest  is  that  found  in  geological  materials  such  as  granular  and  rock,  media,  and 
the  modeling  is  carried  out  using  the  distributed  body  theory  of  Goodman  and  Cowin.  Wave  propagation  is 
studied  through  the  use  of  singular  surface  wave  theory,  and  specific  results  for  the  wave  speed  and  amplitude 
behavior  are  presented.  Three  different  types  of  microstructure  are  modeled  by  using  different  volume 
distribution  functions.  Results  of  wave  speed  and  amplitude  attenuation  are  presented  for  various 
microstructurai  model  parameters. 

1 .  INTRODUCTION 

Geological  materials  such  as  sands,  clays  and/or  rocks  are  complex  materials  and  have  proven  to  be 
difficult  to  model  using  classical  continuum  mechanics.  These  types  of  materials  may  be  classified  as 
materials  with  microstructure  since  at  the  micro-level  the  mass  density  along  with  other  important  field 
variables  are  not  continuous  in  the  mathematical  sense.  Modeling  of  these  materials  using  classical 
continuum  mechanics  (e.g.  elasticity,  plasticity,  viscoelasticity,  etc.)  has  progressed  to  a  point  where 
fundamentally  new  information  will  probably  have  to  come  from  a  theory  incorporating  microstructure  in  its 
basic  framework.  The  work  herein  reported  is  concerned  with  the  modeling  of  the  dynamic  response  of  such 
media.  Of  particular  interest  is  the  propagation  of  mechanical  signals  (i.e.  wave  propagation)  through 
materials  composed  of  granular,  rock,  porous  and  other  discontinuous  structures. 

Studies  of  geological  materials  with  microstructure  started  many  years  ago  with  research  on  granular 
materials  modeled  as  aggregate  assemblies  of  discs  or  spheres.  The  concept  of  modeling  granular  media  as  an 
array  of  elastic  particles  (e.g.  spheres  or  discs)  led  to  the  initial  >  empts  at  predicting  wave  propagation 
through  such  media.  Early  work  [1-6]  employed  a  normal  granular  contact  force  concept  This  initial  work 
investigated  the  propagation  velocity  as  a  function  of  confining  pressure,  particle  size  and  aggregate 
geometry.  It  was  discovered,  however,  that  the  classical  theory  of  contact  due  only  to  normal  forces,  does 
not  in  general  accurately  model  real  materials,  and  thus  Duffy  and  Mindlin  [7]  proposed  a  theory  for  granular 
media  which  included  both  normal  and  tangential  contact  forces.  This  theory  produced  a  non-linear  and 
inelastic  stress-strain  relation. 

More  recent  theories  of  granular  media  behavior  have  included  statistical-stochastic  approaches  [8-13]. 
Cundall  and  Strack  [14]  proposed  a  numerical  method  called  the  distinct  element  technique  for  granular  and 
rock  assemblies,  and  this  approach  has  been  used  for  rubble  screens  [15].  Morland  [16]  considered  a 
rock/granular  media  as  a  regularly  jointed  media  and  used  an  anisotropic  elasticity  approach.  Particulate 
media  has  also  been  studied  by  Hill  and  Harr  based  upon  a  diffusion  equation  derived  from  probabilistic 
models  [17].  Endochronic  theories  have  been  applied  to  granular  soils  [18-20],  and  mixture  theories  [21]  also 
show  some  promise  of  modeling  such  media.  Pore-collapse  models  originally  developed  by  Carrol  and  Holt 
[22]  have  been  used  to  study  the  dynamic  response  of  porous  and  granular  media. 

With  regard  to  experimental  work,  the  method  of  photoelasticity  has  been  employed  to  study  load 
transfer  in  granular  assemblies.  Photoelasticity  has  been  used  to  study  static  behavior  [23-25].  Dynamic 
photomechanics  studies  of  granular  media  have  been  performed  [26-27],  Their  technique  employed  the  use  of 
high  speed  photography  to  record  w^ave  propagation  through  an  assembly  of  birefringent  discs. 
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Some  very  interesting  microstructural  modeling  has  been  done  with  the  so-called  "fabric  tensor"  theories 
for  example  [28].  This  research  has  been  investigating  the  construction  of  constitutive  relations  dependent 
upon  fabric  tensors  which  describe  the  important  microstructural  features  of  the  particular  material,  e.  g„ 
distributions  of  contact  normal  vectors.  At  the  present  time,  wave  propagation  theories  using  fabric  tensors 
have  not  been  developed. 

Finally,  granular  and  porous  materials  have  been  modeled  by  the  so-called  "distributed  body"  theory  [29]. 
This  particular  theory  assumes  that  the  medium  is  distributed  in  space  by  an  independent  kinematical 
function  called  the  volume  distribution  function.  This  theory  has  been  applied  to  wave  propagation  studies 
and  some  success  has  been  achieved  in  modeling  particular  situations.  Based  upon  this  success,  the  present 
work  is  concerned  with  applying  distributed  body  theory  to  wave  propagation  in  geological  soil  media  with 
particular  microstructures. 

The  purpose  of  the  present  work  is  to  extend  the  preliminary  developments  of  the  distributed  body  theory 
for  applications  to  geological  materials.  The  paper  starts  with  the  basic  ideas  previously  developed  for 
one-dimensional  acceleration  waves  propagating  in  a  distributed  body.  Next,  specific  volume  distribution 
functions  are  selected  which  show  promise  of  modeling  granular  geological  materials.  Finally,  a  general 
computer  code  is  developed  which  can  calculate  wave  speed  and  amplitude  behavior  as  a  function  of  distance 
for  all  of  the  various  distribution  functions  selected.  Results  from  the  code  are  shown  to  demonstrate  the 
effects  of  the  selected  volume  distribution  function,  the  constitutive  parameters,  and  the  initial  wave 
amplitude. 

2 .  DISTRIBUTED  BODY  THEORY 


The  distributed  body  theory  originally  developed  by  Goodman  and  Cowin  [29]  was  constructed  to  allow  a 
continuum  theory  to  be  applied  to  materials  with  non-continuous  fields  of  mass  density,  stress,  body  force, 
etc.  Thus,  the  model  could  be  used  to  describe  the  behavior  of  a  wide  variety  of  materials  having  granular 
and/or  porous  structures.  Fundamental  to  the  theory  is  the  assumption  that,  at  any  point  in  the  material,  the 
overall  mass  density  may  be  written  as 

p=vy  (1) 


where  y  is  the  density  of  the  granules  (or  matrix  material)  and  v  =  v(X,t)  is  referred  to  as  the  volume 
distribution  function.  This  function  describes  the  way  thp-medium  is  distributed  in  space  allowing  for  voids 
or  o tiier  particular  granular  structures.  Thus,  this  \hcpky  uncouples  the  mass  density  of  the  granules  from  the 
mass  density  of  the  entire  material,  and  allow^l&^compressibility  due  to  both  granule  compressibility  and 
void  compaction.  In  general  0  <  v  <,  1,  and  v  is  related  to  the  porosity  n  and  void  ^feratio  e  by  the 
expressions 

v  =  1  -  n  =  — !—  .  (2) 


Within  a  one-dimensional  framework,  the  classical  balance  law  of  conservation  of  linear  momentum 


••  9T  .  m 

p  x  =  ^-  +  p  b  (3) 

0  dX  o 

where  T  is  the  stress,  b  is  the  body  force,  x  is  the  particle  position,  X  is  the  reference  position  coordinate, 
and  ( )0  indicate  values  in  the  reference  state.  In  addition  to  this  classical  balance  law,  the  distributed  body 
theory  also  requires  an  independent  balance  equation  governing  the  volume  distribution.  In  one  dimension 
this  second  equation  governing  void  change  is  given  by 

Pok  V  =  H+Pog  (4) 


where  k  is  called  the  equilibrated  inertia,  h  the  equilibrated  stress,  and  g  the  intrinsic  body  force.  Physical 
interpretation  of  the  micro-structural  variables  k,  h  and  g  is  somewhat  difficult  to  make.  In  general,  these 
variables  are  related  to  the  local  contact  mechanics  at  the  granular  level  and  can  be  related  to  particular 
self-equilibrated  singular  stress  states  from  classical  elasticity  (e.  g.  double  force  systems,  centers  of 
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dilatation).  It  has  been  proposed  [31]  that  k  is  related  to  the  void  mean  surface  area  and  to  the  number  of 
voids  present,  h  is  a  result  of  the  interaction  forces  between  neighboring  voids  and  will  vanish  when  the 
voids  are  sufficiently  sepand;  and  g  is  related  to  the  coupling  between  the  total  deformation  of  the  medium 
and  the  changes  in  void  volume. 

For  granular  geological  materials,  we  assume  that  the  media  is  composed  of  compressible  granules  at 
relatively  high  confining  pressures  so  as  to  prevent  material  flow.  For  this  case  an  appropriate  constitutive 
formulation  would  read 


T=T(vo,v’ax,£)  ’  '  (5) 

and  hence  the  stress  depends  upon  the  reference  and  current  volume  distributions,  the  gradient  of  the  volume 
distribution,  and  the  strain  e.  An  explicit  form  of  relation  (5)  which  has  been  proposed  [32]  uses  an  even 
quadratic  form  in  the  gradient  of  v,  i.e. 


T~v  A(vo,v,e)+  1/2 a (vq,  v,e)|  -^-J  •  (6) 

The  constitutive  dependence  on  the  gradient  of  the  volume  distribution  dv/dX  is  significant  and  allows 
an  equilibrium  stress  to  depend  on  this  gradient.  Since  equation  (6)  involves  the  square  of  the  gradient,  it 
will  be  an  isotropic  form  in  that  variable  (required  by  material  frame  indifference)  and,  hence,  the  stress 
response  will  be  independent  of  the  sign  of  the  gradient.  Also,  the  presence  of  the  gradient  term  allows  the 
theory  to  predict  a  generalized  Mohr-Coulomb  failure  criterion. 

First  and  second  order  moduli  defind  by 

dT  1 9v_V 

E=  =v  A  +1/20  I  3X 


▼ 

de2  . 


=?U  +1/2a  |v2 


will  be  needed  for  subsequent  wave  analysis.  Normally  E>0  but  the  second  order  modulus  E,  may  be 
positive  or  negative. 


3.  WAVE  ANALYSIS 


The  basic  premise  of  this  particular  wave  theory  lies  in  modeling  the  wave  as  a  propagating  singular 
surface  across  which  there  exists  a  jump  discontinuity  in  a  particular  variable.  Dynamic  loadings  will 
commonly  produce  either  second-order  acceleration  waves,  having  a  discontinuity  in  the  particle  acceleration 
at  the  wave  front,  or  first-order  shock  wave^  having  a  jump  in  the  particle  velocity  at  the  wave  front. 
Acceleration  waves  will  be  considered  here. 

As  mentioned,  a  wave  is  modeled  as  a  propagating  singular  surface  of  zero  thickness  moving  with  speed 
U.  The  jump  of  a  quantity  <)>  across  this  surface  is  defined  by  [<(>]  =  $'  -  4>+  where  <{>+  and  4>‘  are  the  limiting 
values  of  4>  immediately  ahead  of  and  behind  the  wave,  respectively.  An  acceleration  wave  is  therefore  defined 
as  a  wave  across  which  the  particle  velocity,  strain,  and  volume  distribution  are  continuous  but  their  spatial 
and  temporal  derivatives  are  not.  Thus/this  type  of  motion  carries  propagating  discontinuities  in  the  particle 
acceleration  and  various  other  gradients  of  the  strain  and  volume  distribution.  The  jump  in  the  particle 
acceleration  l[x]  is  called  the  wavq^amplitude,  and  will  be  denoted  by  a(t).  Note  that  for  compressive  waves, 
a(t)  >  0,  whileTor  expansive  waves,  a(t)  <  0. 

Following  singular  surface  analysis  procedures  which  have  now  become  somewhat  standardized,  specific 
relations  for  the  wave  speed  and  amplitude  behavior  can  be  determined.  Specifically,  Nunziato  and  Walsh 
[30]  showed  that  for  distributed  body  theory,  the  wave  speeds  are  given  by  the  roots  of  a  quartic  equation  thus 
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implying  the  existence  of  two  types  of  waves  with  speeds  given  by 


/ 


Up=  1  ^[cJ  +  Cj  +  ^cJ-C*)2  +4(3  j 
l4  =  1/2  [  Cj  +  C\  -J(C]  -  C22)2  +  4p] 

where 


(8) 


(vT)+ 

_ £_ 


(h_)+(T  y 


(9) 


with  subscripts  e,  v  and  X  meaning  partial  differentiation  with  respect  to  the  indicated  variable,  and  ( )+ 
meaning  immediately  ahead  of  the  wave.  The  speed  Up  denotes  the  "fast"  wave  speed  which  is  associated 
predominantly  with  the  elasticity  of  the  granules,  and  will  disappear  if  the  granules  are  incompressible.  The 
quantity  Ug  is  the  "slow"  wave  speed  which  is  connected  to  the  compressibility  of  the  material  due  to 
consolidation,  and  will  vanish  if  there  are  no  voids  (v  =  1). 

With  regard  to  the  wave  amplitude  behavior,  Nunziato  and  Walsh  [30]  have  found  that  the  amplitude  for 
one-dimensional  wave  propagation  satisfies  the  following  nonlinear  Bernoulli  equation 


-^■=  k  (X)  a2  -  n  (X)  a 


(10) 


where  p(X)  and  k(X)  are  material  coefficients  given  in  general  by  rather  lengthy  expressions.  The  coefficient 
p(X)  is  related  to  dispersive  effects,  while  k(X)  reflects  both  the  elastic  response  of  the  granules  and 
dispersive  effects.  Depending  on  the  nature  of  K  and  p,  the  theory  can  predict  growth  or  decay  of  wave 
amplitude.  Using  the  specific  constitutive  form  given  by  Eq.  (6),  for  the  case  of  a  "fast"  wave  propagating 
into  material  at  rest  in  its  reference  condition,  the  coefficients  become 


4(X)  =  -^ 

2Up 


u:  (a) 
— L+ _ L2. 

v  2y 


a v 

ax2 


9v 


ax 


(id 


K(X)  =  . 


-1 


rav\2 

(A  )  +  1/2  (a  ) 

££  0  EE  0 

kaxJ0 

2vyU 


where  the  fast  wave  speed  is  given  by 


V 


with  E0  and  E0  being  the  reference  values  of  the  first  and  second  order  moduli. 
In  order  to  have  real  wave  speeds,  E0  >  0  and  therefore  Eq.  (7)  implies  that 

'id 


(A)  +1/2  (a) 

e  o  e  o 


[dXjo 


>0 


(12) 


(13) 


which  can  be  regarded  as  an  equation  restricting  certain  constitutive  and  micros  true  tural  parameters. 


Journal  of  Wave-Material  Interaction,  Vol.  3,  No.  4,  October  1988 


A4 


Sadd  and  Hossain 


Wave  Propagation  in  Distributed  Bodies  ... 


4.  VOLUME  DISTRIBUTION  FUNCTIONS 

In  order  to  apply  the  distributed  body  theory  and  develop  a  wave  propagation  analysis,  it  is  necessary  to 
have  explicit  constitutive  forms,  see  for  example  Eq.  (6),  and  the  initial  volume  distribution  v0(X)  must  also 
be  specified.  Any  proposed  volume  distribution  function  should  reflect  the  density  variations  and  other 
microstructural  features  within  the  material.  It  is  difficult  to  construct  such  a  function  which  characterizes 
these  variations  precisely  and  yet  has  the  smoothness  requirements  to  be  compatible  with  the  theory.  We 
will  follow  the  approach  that  v0(X)  should  be  a  continuous  function  in  order  to  perform  certain  required 
differentiations  and  integrations  and  'hat  it  yield  the  correct  average  density. 

As  discussed  previously,  in  constructing  a  wave  propagation  analysis,  Nunziato,  et  al.  developed  a 
specific  volume  distribution  function  [32],  Their  work  was  for  a  granular  material,  PBX-9404,  an  explosive 
powder/binder  system.  They  proposed  a  periodic  structure  of  the  form 

V  (X)  =  V  +  (1  —  V  )  COS  — •  (14) 

o  a  *  l 


where  va  and  /  are  material  constants. 

The  quantity  va  would  be  given  by  the  overall  average  density  of  the  material  divided  by  the  granule 
density  and  is  thus  related  to  the  average  value  of  the  volume  distribution.  The  second  material  constant  /  is 
referred  to  as  a  characteristic  length  associated  with  this  periodic  structure.  Clearly  /  specifies  the  length  of 
the  repeating  units  of  the  microstructure.  For  granular  materials,  /  would  be  related,  but  not  necessarily 
equal,  to  the  average  grain  size.  In  regard  to  this  characteristic  length,  the  work  cf  Shahinpoor  [33]  is 
appropriate  to  consider.  Shahinpoor  did  experiments  of  randomly  packed  spherical  granules  on  a  flat  surface. 
His  work  demonstrated  the  concept  of  distinct  packing  geometries  referred  to  as  "Voronoi  cells."  It  is  evident 
that  for  some  packing  geometries,  if  a  periodic  structure  is  assumed,  the  characteristic  length  /  would  be  equal 
to  the  Voronoi  cell  size,  and  thus  could  be  several  grain  diameters. 

Since  the  mechanical  response  of  most  geological  materials  like  sand  or  gravel  is  affected  by  in  situ 
conditions  such  as  overburden,  the  microstructure  will  be  globally  nonhomogeneous,  i.e.,  be  depth 
dependent.  With  this  in  mind,  another  volume  distribution  function  which  can  predict  such  a  structure  may 
be  written  as 


v  (X)=  1  -(1  -v  )e~BX 

0  D 


(15) 


where  vb  and  B  are  material  constants.  Clearly  for  this  case,  the  material  becomes  more  dense  with  depth  X 
into  the  medium.  The  constant  vb  is  the  volume  distribution  at  the  free  surface  X  =  0,  and  the  constant  B 
determines  the  rate  of  consolidation  with  depth.  It  should  be  pointed  out  that  this  exponential  form  does  not 
contain  any  periodic  structure;  hence,  it  should  produce  monotonic  results  for  the  wave  propagation 
characteristics. 

A  final  volume  distribution  function  which  is  proposed  involves  the  combination  of  the  periodic  form 

(14)  and  the  exponential  form  (15).  The  combined  form  employs  simly  the  product  of  these  two  relations, 

A 


v  (X)  = 

o 


„  x  2kX 
Vs  +  (1  -Va)COS-y- 


1-0—  vK)e 


-BX 


(16) 


and  again  va,  vb,  /  and  B  are  material  constants.  It  is  evident  that  this  form  will  thus  produce  a  combined 
periodic -exponential  depth  dependent  microstructure.  The  three  volume  distribution  functions^)  given  by 
Eqs.  (14),  (15)  and  (16)  are  shown  in  Fig.  1. 


4 .  WAVE  MOTION  RESULTS 


Based  upon  the  wave  motion  theory,  a  computer  code  was  developed  to  handle  any  of  the  three  volume 
distribution  functions  given  in  Eqs.  (14),  (15)  and  (16).  The  constitutive  form  employs  Eq.  (6),  with<@ 
specific  values  for  the  two  material  functions  a  and  A  and  their  needed  derivatives.  With  regard  to  the  wave 
speed  calculations,  since  the  variation  of  speed  within  individual  grains  is  normally  not  measurable,  the 
developed  code  calculates  the  average  wave  speed.  This  is  done  through  an  integration  process  over  spatial 
distance  employing  a  numerical  Gauss  quadrature  scheme.  In  order  to  calculate  the  amplitude  behavior,  the 
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non-linear  differential  equation  (10)  was  numerically  integrated  using  a  fourth-order  Runge-Kutta  procedure. 

Specific  wave  motion  results  for  each  of  the  three  proposed  volume  distribution  functions  will  now  be 
presented.  Results  include  wave  speed  and  amplitude  attenuation.  Depending  upon  which  volume 
distribution  function  is  used  in  the  modeling,  a  wide  variety  of  predicted  output  can  result  by  varying 
particular  constitutive,  microstructural,  and  other  input  parameters.  The  constitutive  parameters  in  the  theory 
are  aE,  a^,  Ae  and  Agg,  while  the  microstructural  parameters  are  v,  v^,  Vxx.  va>  vb>  ^  anc*  input 

p  irameter  would  be  the  initial  wave  amplitude  a,,.  In  principle,  the  constitutive  parameters  would  be  a 
function  of  the  microstructural  parameters  as  demonstrated  in  Eq.  (6).  However,  explicit  relationships  for 
these,  parameters  have  not  yet  been  determined.  Therefore,  only  constant  values  will  be  used  for  these  various 
parameters.  /C‘‘  ^ 

For  the  periodic  distribution  model,  the  microstructural  parameters  are  the  average  porqsity  va  and  the 
characteristic  length  /.  Specific  constitutive  parameters  which  /seem  to  give  reasonable  results  for 
geomechanics  materials  are  a£  =  -  450  lb,  Ogg  =(rofFlb,  Ae  =  3  x  (O^db/in2  and  Agg  =  -(£o$  Ib/in. 

Figure  2  illustrates  the  effect  of  porosity  on  the  wave  speed  versus  propagation  depth.  Three  different 
values  of  va  are  shown  using  a  characteristic  length  of  /  =  0.1  in.  As  expected,  the  average  wave  speed 
decreases  with  increasing  porosity.  Figure  3  shows  the  effect  of  the  characteristic  length  on  the  wave  speed, 
and  it  is  seen  that  the  wave  speed  will  increase  with  /.  This  result  is  apparently  related  to  the  fact  that,  with 
an  increase  m  l,  the  wave  will  set  fewer  microstructural  changes  per  unit  length  of  travel  and  will  thus  suffer 
fewer  scatterings.  For  the  periodic  microstructural  model,  the  average  wave  speed  will  oscillate  during  the 
first  few  grains  and  will  then  be  essentially  constant  thereafter. 

Amplitude  behavior  is  illustrated  in  Fig.  4  for  various  cases  of  microstructure  and  initial  input 
acceleration  amplitude.  Figure  4a  shows  the  amplitude  behavior  for  three  different  initial  amplitudes, 
ay  =  5  x  103,  1  x  104,  and  5  x  104  in/s2.  Clearly,  the  expected  result  can  be  seen  in  that  higher  initial 
amplitudes  decay  faster  than  the  lower  amplitude  waves.  Comparing  Figs.  4a  and  4b  indicates  the  effect  of 
porosity  on  the  amplitude  attenuation.  It  is  observed  that  the  attenuation  rate  is  dependent  upon  va  and  as  va 
decreases  (i.e.,  increasing  porosity)  the  rate  of  attenuation  increases.  This  result  is  also  consistent  with  the 
variation  in  wave  speed  with  va  shown  in  Fig.  2.  Comparing  Figs.  4b  and  4c  demonstrates  the  effect  of  /  on 
amplitude  attenuation.  These  fibres  indicate  that  larger  values  of  l  result  in  less  attenuation,  which  is 
consistent  with  the  previous  observjion  regarding  the  variation  of  wave  speed  with  l.  !C  L 

For  the  exponential  volume  distribution  mode)  the  constitutive  parameters  were  chosen  as  otf ■  =  -8  x(£Sr 
Ib,  cijj  =  4.8  x  108  lb,  AE  =  3  x  1(P  lb/in2  and  Agg  =  -1.0  lb/in2.  This  model  contains  the  microstructural 
parameters  of  the  free  surface  porosity  vb  and  the  depth  rate  of  consolidation  B.  As  before,  the  input 
parameter  is  the  initial  amplitude  a^  Figures  5  and  6  show  typical  results  concerning  the  effects  of  these 
parameters  on  the  wave  propagation  variables.  Figure  5  shows  the  effect  of  vb  on  the  average  wave  speed. 
For  this  case,  the  wave  speed  increases  with  depth  due  to  the  overall  decrease  in  porosity  with  depth.  Figure 
6  shows  the  effect  of  vb  and  the  initial  amplitude  on  wave  attenuation.  These  results  give  trends  similar  to 
the  previous  observations  for  the  periodic  distribution  function.^  higher  initial  amplitude  waves 

attenuate  faster  and  the  attenuation  rate  increases  with  porosity.  4  e  - 

For  the  combined  periodic -exponential  distribution  model,  the  chosen  constitutive  parameters  that  were 
used  are  ae  =  73  x  1(H  lb,  aEE  =  5  x  108  lb,  Ae  =  3  x  105  lb/in2  and  AE£  =  -750  lb/in2.  This  model 
contains  aUJdur  microstructural  parameters  va,  vb,  /,  and  B.  These,  along  with  the  initial  amplitude  a,,, 
provide£i6onsiderable  parameter  variations.  Only  a  portion  of  the  possible  parametric  variations  will  te 
presented,  and  these  are  shown  in  Figs.  7  and  8.  Figure  7  shows  the  effect  of  vb  on  the  average  wave  speed 
illustrating  again  a  slower  propagational  speed  for  more  porous  media.  Figure  8  shows  the  effect  of  the 
initial  amplitude  and  porosity  on  wave  attenuation.  These  results  indicate  the  same  trends  as  observed  for  the 
preceding  cases,  namely  an  increase  in  attenuation  with  initial  acceleration  and  porosity. 

5.  SUMMARY  AND  CONCLUSIONS 

This  paper  has  presented  a  wave  propagation  study  based  upon  a  distributed  body  model  of  geological 
media.  A  one-dimensional  theory  has  been  developed  for  explosive  transient  loadings  using  singular  surface 
wave  theory.  General  formulas  for  the  wave  speed  and  amplitude  attenuation  were  taken  from  the  previous 
work  of  Nunziato  ct  al  [30-32],  Specific  forms  for  the  volume  distribution  were  constructed  to  model 
geological  materials,  and  these  were  incorporated  into  the  general  wave  theory.  A  general  computer  code  was 
written  to  calculate  specific  results  from  the  model  for  a  variety  of  constitutive  and  microstructural  model 
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parameters.  Three  specific  forms  of  the  volume  distribution  function  included  a  periodic  model  giving  a 
repeatable  microstructure  with  distance,  an  exponential  form  with  a  depth  dependent  effect  corresponding  to 
ii^feitu  over-burden,  and  a  combined  periodic-exponential  form  yielding  both  repeatable  and  depth  dependent 
microstructure. 

In  regard  to  wave  speed,  since  the  local  speed  through  each  granule  is  difficult  to  observe,  the  average 
wave  speed  was  computed  by  the  code.  It  should  be  pointed  out  that  the  actual  wave  speed  will  vary  in  an 
oscillatory  manner  for  volume  distributions  which  have  a  periodic  nature.  It  is  observed  that  for  the  periodic 
and  periodic -exponential  distribution  cases,  the  average  wave  speed  will  oscillate  through  the  first  couple  of 
periodic  lengths  (barely  discemable  from  Figs.  2,  3,  and  7)  and  then  will  gradually  approach  an  asymptotic 
value.  As  expected^esults  indicate  that  porosity  generally  decreases  the  wave  speed.  However,  it  is  more 
interesting  to  note  that  other  microstructural  parameters  also  have  a  major  effect.  For  example,  results  for 
the  periodic  volume  distribution  indicate  that  the  periodic  length  /,  plays  the  role  of  decreasing  the  wave 
speed  as  l  gets  smaller. 

Although  the  theory  is  capable  of  predicting  both  growth  or  decay  of  wave  amplitude,  for  the  cases 
considered  here  the  wave  amplitude  was  found  to  decrease  with  propagational  distance.  The  attenuation  rate 
was  found  to  be  higher  in  more  porous  media,  and  initially  higher  amplitude  waves  decayed  faster  than  waves 
with  a  lower  initial  amplitude. 

The  work  described  here  deals  solely  with  deterministic  analyses.  Probabilistic  studies  [34]  have  been 
conducted  within  the  framework  of  the  distributed  body  theory.  Such  work  investigated  the  effects  of 
allowing  the  characteristic  length  /,  and  the  average  porosity  va,  for  the  periodic  model,  to  be  random 
variables,  i.e.  to  have  an  average  value  and  a  standard  deviation.  Output  probabilistic  results  were  given  for 
the  expected  value  and  for  plus  or  minus  one  standard  deviation  bounds. 

The  theoretical  results  of  this  work  appear  to  generally  match  both  qualitatively  and  quantitatively  data 
on  wave  propagation  in  soils.  Unfortunately,  experimental  data  on  wave  propagation  through  structured 
geological  media  with  precisely  known  microstructure  was  not  available  to  the  authors.  Consequently, 
comparisons  with  data  could  not  be  made  at  the  present  time.  Current  work  is  under  way  at  the  University  of 
Rhode  Island  to  collect  such  data,  and  future  comparisons  will  be  made. 
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ABSTRACT 

A  computational  and  experimental  study  has  been  conducted  on  the 
propagation  of  mechanical  waves  in  granular  materials.  The  computational 
study  employed  the  use  of  the  distinct  element  method  whereby  the  motion 
of  each  granule  in  the  material  is  modeled  by  rigid-body  dynamics 
assuming  each  particle  interaction  has  a  particular  stiffness  and  damping. 
The  experimental  investigation  has  used  the  method  of  dynamic 
photoelasticity  to  collect  photographic  data  which  provide  information  on 
the  wave  speeds,  inter- granular  contact  loadings,  and  wave  spreading 
characteristics.  The  experimental  results  provide  special  dynamic  material 
constants  necessary  for  the  computational  modeling,  and  they  also  provide 
data  for  comparison  purposes.  Results  from  both  the  computational  and 
experimental  studies  indicate  that  local  microstructure  plays  an  important 
role  in  the  wave  propagation  through  such  materials. 

INTRODUCTION 

A  granular  medium  can  be  characterized  as  a  collection  of  distinct 
particles  which  can  displace  independently  from  one  another  and  interact 
only  through  contact  mechanisms.  Because  of  this  discrete  character,  the 
mechanical  behavior  of  such  materials  under  static  and  dynamic  loading 
conditions  is  very  difficult  to  model.  It  is  now  generally  accepted  that  the 
local  microstructure  or  fabric,  i.e.  the  local  geometrical  arrangement  of 
particles,  plays  a  dominant  role  in  the  transmission  of  mechanical  loadings 
through  these  materials.  Porosity,  which  provides  only  an  average 
estimate  of  microstructure,  is  by  itself  not  sufficient  to  accurately  predict 
the  behavior  of  granular  materials.  Our  aim  here  is  to  understand  the 
dynamic  behavior  of  this  type  of  material  when  it  is  subjected  to 
explosive  loadings  of  short  duration  which  produce  propagating  stress 
waves.  The  discrete  medium  will  act  as  a  structured  wave  guide, 
providing  selective  paths  for  the  waves  to  propagate.  Amplitude 
attenuation  will  then  depend  strongly  upon  the  selected  path  of 


B  1 


propagation,  and  thus  the  wave  propagation  is  linked  to  the  medium 
microstructure. 

A  large  volume  of  reported  research  on  the  mechanical  behavior  of 
granular  materials  exists  in  the  literature.  Constitutive  models  have 
employed  for  example:  elastic/piastic  contact  theories  [1-3],  fabric  tensors 
[4],  distributed  body  models  [5,6],  endochronic  theories  [7],  pore  collapse 
mechanisms  [8],  and  probabilistic  approaches  [9,10].  In  addition,  a 
numerical  scheme  developed  by  Cundall  and  Strack  [11],  called  the 
distinct  element  method  has  also  been  used  to  simulate  granular  media  by 
modeling  the  behavior  of  large  assemblies  of  circular  disks.  In  this 
method,  the  contact  forces  and  displacements  of  an  assembly  of  disks  are 
determined  through  a  series  of  calculations  tracing  the  movements  of  each 
of  the  individual  disks.  The  method  is  based  on  the  use  of  an  explicit 
numerical  scheme  in  which  the  interaction  of  the  granules  is  modeled 
using  rigid-body  dynamics  assuming  each  particle  interaction  has  a 
particular  stiffness  and  damping.  Several  successful  applications  of  this 
method  have  been  reported  [12,13],  and  based  upon  these,  this  method  has 
been  developed  and  applied  to  the  wave  propagation  problems  to  be 
reported  here. 

For  applications  to  wave  propagation,  the  movements  of  each  of  the 
disks  is  a  result  of  the  propagation  through  the  medium  of  disturbances 
originating  at  the  loading  points.  Consequently,  the  wave  speed  and 
amplitude  attenuation  will  be  a  function  of  the  physical  properties  of  the 
discrete  medium,  i.e.  the  microstructure.  Through  the  use  of  base  line 
experimental  data  from  dynamic  photomechanics  studies  on  a  simple 
straight  chain  of  disks,  the  required  dynamic  stiffness  and  damping 
parameters  were  determined.  These  values  were  then  used  to  predict  the 
wave  motion  in  other  more  complex  geometries. 

DISTINCT  ELEMENT  MODELING 

The  distinct  element  approach  uses  Newtonian  rigid-body  mechanics  to 
model  the  translational  and  rotational  motion  of  each  disk  in  a  model 
assembly.  In  the  numerical  routine,  time  steps  are  taken  over  which 
velocities  and  accelerations  are  assumed  to  be  constant.  In  addition,  it  is 
also  assumed  that  during  this  time  step,  disturbances  cannot  propagate 
from  any  disk  further  than  its  immediate  neighbors.  This  then  makes  the 
method  explicit,  and  therefore  at  all  times  the  resultant  forces  on  any  disk 
are  determined  solely  by  its  interactions  with  the  disks  it  is  in  contact. 

Consider  the  case  of  two  disks  in  contact  as  shown  in  Figure  1.  The 
position,  velocity,  acceleration,  angular  velocity,  angular  acceleration, 
radius,  and  mass  of  disk  1  are  labeled  as:  r„  v„  a„  (0„  a„  R„  and  m„ 
with  like  notation  for  disk  2.  The  unit  normal  vector  n  and  unit 
tangential  vector  t  are  defined  as  shown. 


Figure  1.  Schematic  of  Disk  Interaction 

The  normal  component  of  relative  contact  velocity  between  the  two 
disks  is  given  by 


vB  =  (v,  -  v2)-n  (1) 

while  the  tangential  relative  velocity  is 

v,  =  (v,  -  Vj)-t  -  (a>,R,  +  tOjRj)  .  (2) 

Using  a  finite  difference  scheme  with  constant  properties  over  the 
time  interval,  the  relative  velocities  may  be  integrated  with  respect  to  time 
to  yield  the  incremental  relative  normal  and  tangential  displacements,  i.e. 


Axn  =  v„At  =  [(v,  -  Vj)-n]At 
Ax,  =  vAt  =  [(v,  -  Vj)*t  -  (C0[R,  +  Q)2RI)]At  . 


(3) 


In  a  similar  way,  the  absolute  velocity  may  be  computed  from  the 
acceleration  using  the  relation 


Av  =  a  At  .  (4) 

These  relative  displacement  increments  are  to  be  used  with  a 
particular  contact  force-displacement  law  in  order  to  calculate  the  forces 
on  each  disk  in  the  assembly.  Through  allowable  deformations,  the  disks 
in  contact  are  permitted  to  overlap  with  one  another  such  that  the  distance 
between  their  centers  will  become  less  than  (R,  +  Rj).  While  the  general 
technique  could  include  a  complex  nonlinear  contact  law,  the  present  study 
incorporates  a  simple  linear  relation  of  the  form 


AF„  =  KnAxn 
AF,  =  K,Ax,  , 


(5) 
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where  K,  and  K,  are  the  normal  and  tangential  contact  stiffnesses.  At 
each  time  step,  the  force  increments  AFa  and  AF,  are  added  to  the  sum  of 
the  total  forces  Fa  and  F,  on  each  disk  from  previous  time  steps,  i.e. 


(F«)k  =  (FJk,..  +  AFa 
CFJn  =  CFJn-i  +  AF.  , 


(6) 


where  the  indices  N  and  N-l  refer  to  times  tN  and  4,.„  and  At  =  tN-tN.,.  A 
Coulomb-type  friction  law  is  incorporated  to  deal  with  the  tangential 
loading.  This  law  is  defined  by 

(FJ^  =  \i  Fa  +  c  ,  (7) 

where  p.  is  the  coefficient  of  friction  and  c  is  the  cohesion  between  the 
two  disks.  If  the  absolute  value  of  (F,)N  found  from  equation  (6)2  is 
larger  than  (FJ^,  than  (FJ*  is  set  equal  to  (F^. 

The  motion  of  each  disk  is  calculated  using  Newton’s  second  law  of 
motion.  Since  the  behavior  of  real  granular  media  involves  energy 
dissipation,  forms  of  damping  should  be  incorporated  into  the  model. 
Two  forms  of  such  damping  are  therefore  introduced.  A  local  damping 
proportional  to  the  relative  disk  velocities,  and  a  global  damping 
proportional  to  the  absolute  disk  velocities  will  be  included  in  the  force 
balance  laws.  Applying  Newton’s  law  to  disk  1,  therefore  yields 


-  Chvaln  -  C„v„t  -  Cfv,  =  m,a, 
Mi  -  Cttv„R,  -  CgC0]  —  IjOt)  , 


(8) 


where  F  and  M  are  the  resultant  force  and  moment  on  the  disk,  Ctm  and 
C„  are  the  local  damping  coefficients  for  the  normal  and  tangential 
directions,  Cg  is  the  global  damping  coefficient,  and  I,  is  the  moment  of 
inertia  of  the  disk.  Equations  (8)  can  thus  be  solved  for  the  accelerations 
a,  and  a,  over  each  time  increment.  With  the  accelerations  known,  the 
velocities  follow  from  application  of  equation  (4)  and  the  relative 
displacements  can  then  be  computed  from  equation  (3).  This  leads  to  new 
values  of  the  forces  through  (5)  for  the  next  time  increment,  and  the  cycle 
is  repeated  again  for  each  disk.  In  this  manner,  large  assemblies  of  disks 
can  be  analyzed  in  a  reasonable  amount  of  computer  time.  Values  of  the 
stiffness  and  damping  parameters  appropriate  for  a  given  material  are 
difficult  to  measure.  For  the  static  case,  the  stiffness  properties  may  be 
computed  from  Hertz  theory  or  from  other  elasticity  analyses.  However, 
for  the  dynamic  case  involving  loadings  of  short  duration,  estimates  of  the 
model  parameters  are  extremely  difficult  to  make.  It  is  here  that 
experimental  dynamic  photoelasticity  can  be  used. 
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EXPERIMENTAL  TECHNIQUE 


The  optical  technique  of  dynamic  photoelasticity  along  with  high  speed 
photography  is  very  well  suited  to  study  wave  propagation  in  different 
assemblies  of  granular  media.  Details  on  the  wave  speed,  inter-granular 
contact  forces,  and  wave  spreading  phenomena  may  all  be  determined 
through  this  experimental  method.  The  granular  medium  was  simulated 
by  assemblies  of  one  inch  diameter  discs  which  were  fabricated  from  a 
brittle  polyester  material,  Homalite  100.  The  dynamic  loading  was 
achieved  by  exploding  a  small  charge  of  PETN  in  a  specially  designed 
charge  holder  directly  on  top  of  one  of  the  grains.  Experimental  granular 
assemblies  were  placed  in  the  optical  bench  of  a  high  speed  multiple  spark 
gap  camera.  The  details  of  this  camera  are  given  in  [14].  This  high 
speed  photographic  system  operates  as  a  series  of  high  intensity,  extremely 
short  duration  pulses  of  light  and  provides  20  photoelastic  images  at 
discrete  times  during  the  dynamic  event.  Framing  rates  of  up  to  106 
frames  per  second  are  attainable  using  this  photographic  stem.  Typical 
photographs  obtained  during  the  experiments  are  shown  in  Figures  2  and 
3.  These  photographs  show  the  isochromatic  fringes  at  different  times  as 
the  stress  wave  passes  through  the  granular  assemblies.  Figure  2 
illustrates  the  propagation  of  a  stress  wave  down  a  single  straight  chain  of 
disks,  while  Figure  3  shows  the  dynamic  fringe  patterns  associated  with 
waves  moving  through  three  different  two-dimensional  assemblies. 
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Figure  2.  Photoelastic  Fringe  Patterns  of  the  Single  Disk  Chain 
for  a  Sequence  of  Times 
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FRAME  12,  t  =  125^3 


k 

(b)  Hexagonal  Close  Packing 


FRAME  18.  I  =159.5/i» 

(c)  Random  Packing 
Figure  3.  Photoelastic  Fringe  Patterns 


The  isochromatic  fringes  obtained  during  the  experiments  are  lines  of 
constant  maximum  shear  stress.  These  fringes  are  related  to  the  stress 
field  by  the  stress  optic  law 


a,  -  a2  =  Nfyh  ,  (9) 

where,  and  ct2  are  the  principal  stresses,  N  is  the  fringe  order,  f„  is  the 
material  fringe  value,  and  h  is  the  thickness  of  the  model. 

Shukla  and  Nigam  [15]  have  shown  that  equation  (9)  when 
combined  with  the  Hertz  contact  stress  equations  can  be  used  to  analyze 
the  stress  field  in  the  vicinity  of  the  contacts  shown  in  the  Figures  2  and 
3.  This  analysis  yields  both  the  normal  and  the  tangential  contact  forces 
as  a  function  of  time.  Some  particular  values  of  the  peak  contact  forces 
obtained  in  these  geometries  are  shown  in  Figures  4  and  5,  and  these  are 
to  be  used  in  conjunction  with  the  numerical  study. 

RESULTS  AND  CONCLUSIONS 

In  order  to  apply  the  distinct  element  method  to  various  granular 
assemblies,  values  for  the  model  parameters  (stiffness  and  damping 
coefficients)  must  be  determined.  For  the  dynamic  case,  we  expect  that 
the  contact  stiffnesses  and  K,,  will  be  difficult  to  calculate  from  simple 
static  Hertz  contact  theory.  In  addition,  values  for  the  damping  would 
also  be  difficult  to  calculate  from  say  elastodynamic  theory.  Hence,  the 
dynamic  photoelastic  experiments  were  used  to  provide  this  information. 
Experiments  were  performed  on  a  single  straight  chain  of  disks  (see 
Figure  2),  and  these  experiments  were  considered  to  be  our  calibration 
tests.  Appropriate  values  for  these  stiffness  and  damping  coefficients  were 
thus  determined  to  match  the  data  from  the  calibration  experiments. 
These  values  were  then  retained  in  the  model  to  predict  the  wave  motion 
in  other  more  complicated  geometries.  For  the  cases  to  be  reported  here 
the  contact  stiffnesses  were  taken  as  K,  =  K,  =  6.4  x  106  N/m,  and  the 
local  normal  damping  coefficient  was  =  32  N-s/m.  Global  damping 
was  not  included,  and  tangential  contact  loading  was  set  to  zero.  The 
input  loading  was  chosen  to  be  a  triangular  profile  of  60|is  duration,  and 
the  time  step  was  taken  as  At  =  2ps. 

Results  from  the  distinct  element  modeling  are  shown  in  Figures  4 
and  5  for  four  different  assemblies,  i.e.  the  single  straight  chain,  body- 
centered  cubic,  hexagonal  close  packing,  and  a  random  packing. 
Numerical  predictions  of  the  maximum  inter-granular  contact  loadings 
(normalized  with  respect  to  the  input  loading)  are  shown  at  various 
contacts  in  the  assemblies.  For  comparison,  the  corresponding 
experimental  values  are  also  given  in  parentheses.  As  menu  ,ned,  the 
straight  chain  experiment  was  actually  used  to  determine  the  model 
parameters.  For  the  body  centered  cubic  geometry  shown  in  Figure  4, 
experiments  indicate  that  for  input  excitation  directed  along  a  column 
chain,  the  resulting  wave  motion  will  occur  only  along  that  chain. 
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Single  Disk  Chain  Body-Centered  Cubic  Geometry 


Figure  4.  Distinct  Element  Predictions  and  Comparisons 

This  was  also  predicted  numerically  as  seen  from  Figure  4.  On  the  other 
hand,  the  hexagonal  close  packing  geometry  allows  wave  motion  in  many 
different  chains  and  paths  (this  can  be  seen  from  the  experiments  shown 
in  Figure  3).  Some  of  the  numerical  predictions  are  compared  with 
experiments  for  this  geometry  for  two  principal  chains  as  shown  in  Figure 
5.  The  random  model  material  offers  the  most  complex  wave 
transmission  phenomena.  This  model  contains  some  grains  with  noticeable 
tangential  inter-granular  contact  loadings.  However,  as  seen  in  Figure  5, 
our  frictionless  numerical  simulation  provides  contact  values  which  are 
generally  comparable  to  the  experimental  data. 

Results  of  this  study  illustrate  the  dependence  of  microstructure  on 
the  wave  propagation.  In  order  to  describe  local  microstructure  or  fabric, 
it  is  convenient  to  define  the  term  branch  vectors  which  are  drawn  from 
the  mass  centers  of  adjacent  disks.  The  angle  between  neighboring  branch 
vectors  may  be  denoted  as  a  branch  angle.  It  is  therefore  clear  from 
Figures  4  and  5  that  the  inter-granular  contact  force,  which  is  a  measure 
of  the  wave  amplitude,  behaves  quite  differently  when  the  wave 
propagates  through  media  with  large  changes  in  the  branch  vector 
distribution.  Further  numerical  and  experimental  studies  are  underway  in 
an  attempt  to  model  additional  dynamic  phenomena  in  granular  media. 
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Abstract 

Theoretical  and  experimental  studies  have  been  conducted  on  the  dynamic  response 
of  granular  materials  containing  local  discontinuities  of  voids  and  inhomogeneous 
inclusions.  The  granular  medium  was  simulated  by  a  specific  assembly  of  circular  disks 
which  were  subjected  to  explosive  loadings  of  short  duration.  Voids  were  created  by 
removing  particular  disks  from  the  assembly,  while  inclusions  were  constructed  by 
replacing  certain  disks  with  ones  of  a  higher  impedance  material.  The  computational 
simulation  was  accomplished  through  the  use  of  the  distinct  element  method  in  which  the 
intergranular  contact  forces  and  displacements  of  the  assembly  disks  are  determined  through 
a  series  of  calculations  tracing  the  movements  of  each  of  the  individual  disks.  The 
experimental  study  employed  the  use  of  photoelasticity  in  conjunction  with  high  speed 
photography  to  collect  photographic  data  of  the  propagation  of  waves  in  transparent 
assemblies  of  model  granular  media.  Comparisons  were  made  between  the  computational 
results  and  the  experimental  data  for  the  local  intergranular  contact  forces  around  each  void 
or  inclusion.  Both  voids  and  inclusions  produce  local  wave  scattering  through  various 
reflection  mechanisms,  and  the  results  seem  to  indicate  that  the  inclusions  produced  higher 
local  wave  attenuation. 

1.  Introduction 

Granular  media  can  be  described  as  a  collection  of  distinct  particles  which  can 
displace  independently  from  one  another  and  which  interact  only  through  contact 
mechanisms.  This  type  of  media  transmits  loadings  through  discrete  paths,  and  therefore 
the  mechanical  behavior  of  such  materials  under  static  and  dynamic  loading  conditions  is 
very  complex  and  difficult  to  model.  It  is  now  generally  accepted  that  the  local 
microstructure  or  fabric,  i.e.  the  local  geometrical  arrangement  of  particles,  plays  a 
dominant  role  in  the  transmission  of  mechanical  loadings  through  these  materials.  The 
material  porosity  provides  only  an  average  estimate  of  microstructure,  and  by  itself  it  is  not 
sufficient  to  accurately  predict  the  behavior  of  granular  materials.  Our  general  aim  here  is 
to  understand  the  dynamic  behavior  of  this  type  of  material  when  it  is  subjected  to 
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explosive  loadings  of  short  duration  which  produce  propagating  stress  waves.  The  discrete 
medium  will  act  as  a  structured  wave  guide,,  providing  selective  paths  for  the  waves  to 
propagate.  Amplitude  attenuation  will  then  depend  strongly  upon  the  selected  path  of 
propagation,  and  thus  the  wave  propagation  is  strongly  linked  to  the  medium 
microstructure. 

A  considerable  amount  of  research  on  the  mechanical  behavior  of  granular  materials 
has  been  reported  in  the  literature.  Constitutive  models  have  employed  for  example: 
elastidplastic  contact  theories  (Deresiewicz,  1959,  Walton,  1987,  and  Petrakis  and  Dobry, 
1986),  fabric  tensors  (Nemat-Nasser,  1983),  distributed  body  models  (Goodman  and  Cowin, 
1972),  endochronic  theories  (Bazant,  1983),  pore  collapse  mechanisms  (Carroll  and  Holt, 
1972),  and  probabilistic  approaches  (Endley  and  Peyrot,  1977,  and  Fu,  1977).  The  concept 
of  modeling  granular  media  as  an  array  of  elastic  disks  or  spheres  lead  to  the  initial 
attempts  at  predicting  wave  propagation  phenomena.  Early  wave  propagation  studies 
include  Iida  (1939),  Takahashi  and  Sato  (1949),  Hughes  and  Cross  (1951),  Hughes  and 
Kelly  (1952),  Gassman  (1951),  Brandt  (1955),  and  Duffy  and  Mindlin  (1957).  This  initial 
work  investigated  the  propagation  velocity  as  a  function  of  confining  pressure,  panicle  size 
and  geometrical  packing.  This  particular  modeling  concept  lead  to  work  in  determining  the 
elastic  constants  of  particular  granular  assemblies,  see  for  example,  Hendron  (1963), 
Petrakis  and  Dobry  (1986),  and  Walton  (1987).  In  addition,  wave  propagation  studies  for 
granular  and  porous  media  by  Nunziato  et.al.  (1978)  and  Sadd  and  Hossain  (1989)  have 
employed  the  distributed  body  theory.  Experimental  studies  of  this  problem  employing 
the  method  of  dynamic  photoelasticity  have  been  reported  by  Shukla  et.al. 
(1985,1987,1988).  In  this  method,  high  speed  photography  was  used  to  collect 
photographic  data  of  wave  propagation  in  transparent  model  materials  composed  of 
assemblies  of  birefringent  disks. 

A  numerical  scheme  developed  by  Cundall  and  Strack  (1979),  called  the  distinct 
element  method  has  also  been  used  to  simulate  granular  media  by  modeling  the  behavior  of 
large  assemblies  of  circular  disks.  In  this  method,  the  contact  forces  and  displacements  of 
an  assembly  of  disks  are  determined  through  a  series  of  calculations  tracing  the  movements 
of  each  of  the  individual  disks.  The  metlnd  is  based  on  the  use  of  an  explicit  numerical 
scheme  in  which  the  interaction  of  the  granules  is  modeled  using  rigid-body  dynamics 
assuming  each  particle  interaction  has  a  particular  stiffness  and  damping.  For  applications 
to  wave  propagation,  the  movements  of  each  of  the  disks  is  a  result  of  the  propagation 
through  the  medium  of  disturbances  originating  at  the  loading  points.  Consequently,  the 
wave  speed  and  amplitude  attenuation  (intergranular  contact  force)  will  be  a  function  of  the 
physical  properties  of  the  discrete  medium,  i.e.  the  microstructure.  Several  successful 
applications  of  this  method  have  been  reported  (Thorton  and  Barnes,  1986,  Bathurst  and 
Rothenburg,  1988,  and  Sadd  et.al.,  1989),  and  based  upon  these,  this  method  has  been 
developed  and  applied  to  the  wave  propagation  problems  to  be  reported  here. 

The  present  study  focuses  on  a  specific  aspect  of  wave  propagation  in  granular 
materials,  namely  the  effects  of  voids  and  inclusions.  It  is  well  known  that  actual  granular 
media  contains  both  voids  and  heterogeneous  inclusions.  These  quantities  further 
complicate  an  already  complex  microstructural  material.  Most  past  studies  on  wave 
propagation  in  these  materials  have  been  limited  to  looking  at  aggregate  assemblies  with 
uniform  packing  geometries.  Our  focus  here  is  to  investigate  the  local  effects  produced  by 
voids  and  inclusions  in  regular  arrays  of  circular  disks,  and  of  primary  concern,  is  the 
wave  scattering  in  the  vicinity  of  the  microstructural  defect.  The  wave  propagation 


phenomena  is  to  be  studied  through  the  determination  of  the  intergranular  contact  force 
distribution  in  the  neighborhood  of  a  particular  void  or  inclusion.  This  problem  has  been 
studied  using  the  computational  method  of  distinct  elements,  and  using  the  experimental 
technique  of  dynamic  photoelasticity. 

2.  Distinct  Element  Method 

The  distinct  element  method  is  a  simplified  modeling  concept  which  uses  Newtonian 
rigid-body  mechanics  to  model  the  translational  and  rotational  motion  of  each  disk  in  a 
model  assembly.  The  technique  establishes  a  discretized  time  stepping  numerical  routine, 
in  which  granule  velocities  and  accelerations  are  assumed  to  be  constant  over  each  time 
interval.  It  is  also  assumed  that  during  each  time  step,  disturbances  cannot  propagate  from 
any  disk  further  than  its  immediate  neighbors.  Under  these  assumptions,  the  method 
becomes  explicit,  and  therefore  at  any  time  increment  the  resultant  forces  on  any  disk  are 
determined  solely  by  its  interactions  with  the  disks  it  is  in  contact 

In  order  to  describe  the  method,  consider  the  case  of  two  disks  in  contact  as  shown 
in  Figure  1.  The  position,  velocity,  acceleration,  angular  velocity,  angular  acceleration, 
radius,  and  mass  of  disk  1  are  labeled  as:  r,,  v„  C0i,  a„  R„  and  mlf  with  like  notation 
for  disk  2.  The  unit  normal  vector  n  and  unit  tangential  vector  t  are  defined  as  shown. 

The  normal  component  of  relative  contact  velocity  between  the  two  disks  is  given 
by 


v„  =  (vt  -  v2)-n  (2.1) 

while  the  tangential  relative  velocity  is 

vt  =  (v,  -  VjM  -  (cDjR,  +  ct^Rj)  .  (2.2) 

Using  a  finite  difference  scheme  with  constant  properties  over  the  time  interval,  the 
relative  velocities  may  be  integrated  with  respect  to  time  to  yield  the  incremental  relative 
normal  and  tangential  displacements,  i.e. 

Axn  =  v„At  =  [(v,  -  v2)-n]At 

Ax,  =  v,At  =  [(vt  -  VjH  -  (co,R,  +  (OjRjXIAt  .  (2.3) 

In  a  similar  way,  the  absolute  velocity  may  be  computed  from  the  acceleration  using  the 
relation 


Av  =  a  At  .  (2.4) 

These  relative  displacement  increments  are  to  be  used  with  a  particular  contact 
force-displacement  law  in  order  to  calculate  the  forces  on  each  disk  in  the  assembly. 
Through  allowable  deformations,  the  disks  in  contact  are  permitted  to  overlap  with  one 
another  such  that  the  distance  between  their  centers  will  become  less  than  (R,  +  R2). 
While  the  general  technique  could  include  a  complex  nonlinear  contact  law,  the  present 
study  incorporates  a  simple  linear  relation  of  the  form 
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AFa  =  K^Ax, 
AF,  =  K^Ax,  , 


(2.5) 


where  and  K,  are  the  normal  and  tangential  contact  stiffnesses.  At  each  time  step,  the 
force  increments  AFn  and  AF,  are  added  to  the  sum  of  the  total  forces  Fn  and  F,  on  each 
disk  from  previous  time  steps,  i.e. 


(F„)n  =  (F.V,  +  AFn 

(FDn  =  (F^.,  +  AF,  ,  (2.6) 

where  the  indices  N  and  N-l  refer  to  times  tN  and  tN.,,  and  At  =  tN-tN_,.  A  Coulomb-type 
friction  law  is  incorporated  to  deal  with  the  tangential  loading.  This  law  is  defined  by 

(Fdw  =  4  Fn  +  c  ,  (2.7) 

where  |i  is  the  coefficient  of  friction  and  c  is  the  cohesion  between  the  two  disks.  If  the 
absolute  value  of  (Ft)N  found  from  equation  (2.6)2  is  larger  than  (FJ^,  then  (F,)N  is  set 
equal  to  (FJ^,. 

Using  Newton’s  second  law  of  motion,  the  acceleration  of  each  disk  at  each  time 
interval  can  be  determined.  Now  since  the  behavior  of  real  granular  media  involves  energy 
dissipation,  the  modeling  introduces  damping  mechanisms.  Two  forms  of  such  damping 
are  therefore  introduced.  A  local  damping  proportional  to  the  relative  disk  velocities,  and  a 
global  damping  proportional  to  the  absolute  disk  velocities  will  be  included  in  the  force 
balance  laws.  Applying  Newton’s  law  to  disk  1,  therefore  yields 


F,  -  Q,valn  -  C„v„t  -  C,v,  =  m.a, 

Mi  -  C(lv„R,  -  0,0),  =  I,a,  ,  (2.8) 

where  F  and  M  are  the  resultant  force  and  moment  on  the  disk,  C,„  and  C„  are  the  local 
damping  coefficients  for  the  normal  and  tangential  directions,  C,  is  the  global  damping 
coefficient,  and  I,  is  the  moment  of  inertia  of  the  disk.  Equations  (2.8)  can  thus  be  solved 
for  the  accelerations  a,  and  a,  over  each  time  increment.  With  the  accelerations  known, 
the  velocities  follow  from  application  of  equation  (2.4)  and  the  relative  displacements  can 
then  be  computed  from  equation  (2.3).  This  leads  to  new  values  of  the  forces  through 
(2.5)  for  the  next  time  increment,  and  the  cycle  is  repeated  again  for  each  disk.  In  this 
manner,  large  assemblies  of  disks  can  be  analyzed  in  a  reasonable  amount  of  computer 
time.  Values  of  the  stiffness  and  damping  parameters  appropriate  for  a  given  material  are 
difficult  to  measure.  For  the  static  case,  the  stiffness  properties  may  be  computed  from 
Hertz  theory  or  from  other  elasticity  analyses.  However,  for  the  dynamic  case  involving 
loadings  of  short  duration,  estimates  of  the  model  parameters  are  difficult  to  make.  In 
order  to  determine  estimates  of  these  model  parameters,  experimental  results  from  dynamic 
photoelasticity  were  employed.  Data  from  simple  calibration  tests  on  a  single  chain  of 
disks  were  used  to  determine  values  for  the  contact  stiffness  and  damping  for  the  disk 
assemblies  under  study.  For  the  stiffness  parameter,  the  dynamic  stiffness  was  written  as 

K,  =  a  K,(,), 
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(2.9) 


(2.20) 


where  K^(,>  is  the  static  contact  stiffness  from  Hertz  theory  given  by 

^  (■) _  TthEjE^ 

^  "  2(E,+E2) 

and  a  is  a  dynamic  stiffness  coefficient  which  is  to  be  determined  by  the  experiments,  h  is 
the  disk  thickness,  and  Et  and  Ej  are  the  elastic  modulli  of  the  two  disks  in  contact.  The 
a  coefficient  may  be  thought  of  as  an  adjustment  parameter,  less  than  unity,  which 
accounts  for  the  fact  that  the  entire  disk  will  not  deform  during  the  dynamic  event. 

3.  Photomechanics  Studies 

A  series  of  dynamic  photoelastic  experiments  were  conducted  to  provide 
experimental  information  on  the  effect  of  voids  and  inclusions  on  the  wave  propagation 
phenomena.  For  the  experimental  study,  the  granular  medium  was  simulated  with 
assemblies  of  polymeric  birefringent  disks  of  Homalite  100,  25.4mm  in  diameter  and 
6.25mm  thick.  In  all  the  experiments  reported  here,  the  disks  were  assembled  in  an 
hexagonal  close  packing  geometry  as  shown  in  Figure  2.  Voids  in  the  assembly  were 
created  by  removing  disks  from  different  locations,  while  inclusions  were  created  by 
replacing  particular  disks  with  ones  of  a  different  material  (steel).  Dynamic  loading  was 
achieved  by  detonating  a  small  charge  of  PETN  in  a  specially  designed  charge  holder, 
which  was  mounted  at  the  top-centerline  of  the  model  assemblies.  The  experimental 
models  were  placed  in  the  optical  bench  of  a  high  speed  photographic  system.  This  high 
speed  photographic  system  operated  as  a  series  of  high  intensity,  extremely  short  duration 
pulses  of  light  and  provided  20  photoelastic  images  at  discrete  times  during  the  dynamic 
event.  Framing  rates  of  up  to  10s  frames  per  second  are  achievable,  and  this  allows  studies 
of  wave  propagation  to  be  made,  see  Riley  and  Dally  (1969). 

A  typical  sequence  of  four  images  from  each  experiment  are  shown  in  Figures  2 
through  6.  The  photographic  data  shows  the  isochromatic  fringe  patterns  at  different  times 
as  the  stress  wave  propagates  through  the  model  assemblies.  The  wave  propagation  history 
can  thus  be  clearly  seen  from  a  sequence  of  such  photographs.  Figure  2  illustrates  the 
wave  patterns  for  the  assembly  with  no  voids  or  inclusions.  Figures  3  and  4  show  the 
wave  motion  for  the  cases  of  single  and  multiple  voids,  while  Figures  5  and  6  illustrate  the 
cases  with  inclusions  of  higher  density  and  stiffness.  Inspection  of  the  photographs  reveals 
that  the  wave  length  X,  of  the  loading  pulse  is  much  larger  than  the  disc  diameter  D;  in 
fact  X=4D.  Furthermore,  in  most  cases  the  fringe  pattern  around  the  contact  points  are 
symmetric  on  either  side  of  the  contact  points  and  are  similar  to  the  fringes  obtained  under 
static  compression.  Both  these  features  indicated  that  around  the  contact  zone,  quasi  static 
loading  was  present  during  the  wave  propagation  event.  Thus  Hertz  theory  can  be  used  to 
estimate  the  stress  field  in  the  vicinity  of  the  contacts. 

The  isochromatic  fringes  photographed  during  the  experiments  are  lines  of  constant 
maximum  shear  stress,  and  are  related  to  the  stress  field  by  the  stress  optic  law 

a,  -  o2  =  MJh  (3.1) 

where  a,  and  o2  are  the  principal  stresses,  N  is  the  fringe  order,  f0  is  the  material  fringe 
value,  and  h  is  the  model  thickness.  Using  relation  (3.1)  along  with  the  Hertz  contact 
stress  equations,  Shukla  and  Nigam  (1985)  have  developed  a  scheme  to  compute  the  quasi- 
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static  stress  field  in  the  vicinity  of  the  contacts  between  each  grain.  These  contact  stresses 
may  then  be  integrated  along  the  contact  length  to  determine  the  normal  and  tangential 
contact  forces  transmitted  between  the  disks  at  various  times  during  the  dynamic  process. 
Thus  the  intergranular  contact  loading,  which  is  related  to  the  wave  amplitude,  can  be 
experimentally  computed.  These  experimental  results  can  then  be  compared  with  the 
theory,  and  they  can  also  be  used  to  determine  necessary  material  model  parameters  (e.g. 
the  contact  stiffness  and  damping)  in  order  to  use  the  distinct  element  method. 

4.  Results  and  Comparisons 

Specific  results  of  the  distinct  element  modeling  along  with  comparisons  to  the 
experimental  data  will  now  be  given  for  particular  granular  assemblies  containing  voids  and 
inclusions.  The  basic  assembly  geometry  is  a  hexagonal  close  packing  arrangement  shown 
in  Figure  7.  The  theoretical  model  used  an  input  loading  of  triangular  time  dependence 
with  a  60p.s  duration  to  model  the  explosive  loading  from  the  experiments.  The  local 
stiffness  and  damping  parameters  K^,  K,,  C,„,  and  Cu  needed  in  the  distinct  element  model 
were  determined  from  experimental  tests  conducted  on  straight  single  chains  of  disks. 
These  calibration  experiments  provided  specific  values  for  the  wave  speed  and  amplitude 
attenuation  for  a  very  simple  geometry.  The  model  parameters  were  thus  chosen  to 
produce  the  best  match  to  the  calibration  data,  and  these  values  were  then  retained  in  the 
model  to  be  used  for  calculations  of  more  complicated  geometry.  The  contact  stiffnesses 
for  Homalite-to-Homalite  contacts  were  taken  as  =  K,  =  6.4  x  106  N/m  (a  =  0.27),  and 
the  local  normal  damping  coefficient  was  Cta  =  32  Ns/m.  Global  damping  was  not 
included  ,  and  tangential  contact  loading  was  set  to  zero.  The  time  stepping  increment  was 
taken  to  be  2p.s  for  all  cases  studied,  and  this  was  felt  to  be  appropriate  to  calculate  the 
essential  features  of  the  dynamic  event 

Results  from  the  distinct  element  technique  are  shown  in  Figure  7  for  the  basic 
granular  media  with  no  voids  or  inclusions.  Numerical  predictions  of  the  maximum  inter¬ 
granular  contact  loadings  (normalized  with  respect  to  the  input  loading)  are  shown  at 
various  contacts  in  the  assembly.  The  corresponding  experimental  values  are  also  given  in 
parentheses.  The  contact  loading  values  are  symmetric  about  the  assemblies’  vertical 
centerline,  and  the  results  indicate  the  rapid  attenuation  which  occurs  along  various  chains 
or  paths  in  the  assembly.  Comparison  of  theoretical  with  experimental  contact  loadings 
indicate  that  they  differ  by  an  average  amount  of  25%.  Experimental  determination  of  the 
wave  speed  was  accomplished  from  the  known  position  of  the  fringe  patterns  in  the 
photographs.  These  results  indicated  that  the  leading  wave  front  in  the  Homalite 
assemblies,  propagates  at  approximately  995  m/s,  which  is  about  50%  of  the  P-wave  speed 
in  the  virgin  disk  material.  Wave  speed  predictions  from  the  distinct  element  model 
matched  well  (within  10%)  to  these  measured  values. 

Figures  8  and  9  illustrate  the  case  of  a  granular  media  with  voids  present.  Figure  8 
contains  a  single  void,  while  Figure  9  contains  a  series  of  voids  down  the  vertical 
centerline  of  the  assembly.  The  presence  of  a  void  produces  significant  local  wave 
scattering  and  attenuation,  especially  along  the  main  vertical  chain  down  the  centerline  of 
the  assembly.  Local  inter-granular  contact  forces  become  elevated  near  the  void;  however, 
at  contacts  remote  from  the  void,  the  loading  values  are  similar  in  magnitude  (with  the 
exception  of  points  on  the  vertical  centerline)  to  those  in  Figure  7.  Computational  and 
experimental  contact  loads  compare  to  within  an  average  difference  of  15%  for  the 
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assemblies  with  voids. 

The  cases  of  granular  media  with  inclusions  of  different  material  are  shown  in 
Figures  10  and  11.  This  situation  is  attempting  to  model  the  local  effects  of  heterogeneity 
in  a  granular  medium.  Figure  10  illustrates  the  case  of  a  single  inclusion,  and  Figure  1 1 
contains  the  case  of  several  inclusions  down  the  vertical  centerline  of  the  assembly.  As 
mentioned  the  inclusion  is  simply  a  disk  of  a  different  material,  and  in  this  case  the 
inclusion  material  was  steel,  which  in  comparison  to  the  granular  material  (Homalite  100) 
has  a  much  higher  stiffness  and  density.  The  contact  stiffness  between  the  Homalite  and 
steel  disks  was  calculated  using  relations  (2.9)  and  (2.10)  retaining  the  same  value  of  a. 
The  local  damping  parameter  C,„,  was  kept  the  same  for  this  case,  since  there  was  only  a 
small  number  of  inclusions  present. 

There  will  be  an  impedance  mismatch  at  the  contacts  of  the  Homalite  and  steel,  and 
thus  there  will  be  a  sizeable  difference  between  the  reflection  and  transmission  phenomena 
at  these  contacts  in  comparison  to  those  of  the  rest  of  the  medium.  For  example,  at  point 
C  in  Figure  10,  the  wave  is  attempting  to  propagate  from  a  relatively  soft  material  into  a 
stiff  material.  Consequently,  a  sizeable  reflection  occurs  at  this  contact  producing  an 
upward  traveling  wave  and  a  large  contact  force.  Very  little  wave  motion  is  transmitted 
into  the  inclusion,  and  thus  the  inclusion  acts  to  block  the  wave  motion  along  the  vertical 
disk  chain.  For  the  inclusion  cases,  the  average  difference  between  theoretical  predictions 
and  experimental  data  was  20-30%. 

5.  Conclusions 

The  methods  of  distinct  elements  and  dynamic  photomechanics  have  been  used  to 
study  the  effects  of  voids  and  heterogeneous  inclusions  on  the  wave  propagation  in 

granular  materials.  The  granular  medium  was  simulated  by  a  specific  assembly  of  circular 

disks  arranged  in  an  hexagonal  close  packing  geometry.  The  voids  were  created  by 

removing  particular  disks,  while  inclusions  were  constructed  by  replacing  particular  disks 
with  ones  of  a  higher  modulus  material.  Comparisons  were  made  between  the 

computational  results  and  the  experimental  data  for  the  local  intergranular  contact  forces 
around  each  void  or  inclusion.  Although  comparisons  produced  some  differences  as  high 
as  30%,  it  is  felt  that  since  the  experimental  data  itself  contains  scatter  of  approximately 
10%,  the  computational  scheme  does  provide  reasonable  predictions.  Improvements  of  the 
modeling  procedures  could  be  accomplished  by  incorporating  a  more  sophisticated  dynamic 
contact  law. 

Comparing  the  results,  it  is  apparent  that  both  the  voids  and  inclusions  cause  local 
wave  scattering.  In  comparison  to  the  case  with  no  voids  or  inclusions,  elevated  contact 
forces  occur  near  the  discontinuity,  and  along  particular  paths  dictated  by  the  local 
microstructure,  rapid  wave  amplitude  attenuation  occurs.  A  void  produces  wave  scattering 
through  ffee-surface  reflection  from  the  empty  volume,  whereas  the  inclusion  causes 
sizeable  reflections  from  the  material  of  higher  impedance.  Comparing  the  path  .  AB  in 
Figures  8  and  10,  it  appears  that  the  inclusion  produces  higher  attenuation  than  the  void. 
Current  work  continues  on  these  studies,  and  investigations  are  underway  on  wave 
propagation  in  granular  materials  with  additional  and  more  complex  microstructural 
features. 
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Figure  1.  Schematic  of  Disk  Interaction 
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Figure  7.  Normalized  Maximum  Contact  Loading  Comparisons:  No  Voids  or  Inclusions. 


Experimental  Data  in  Parentheses 


C15 


Figure  8.  Normalized  Maximum  Contact  Loading  Comparisons:  Single  Void. 
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Figure  9.  Normalized  Maximum  Contact  Loading  Comparisons:  Series  of  Voids. 
Experimental  Data  in  Parentheses 


Figure  10.  Normalized  Maximum  Contact  Loading  Comparisons:  Single  Inclusion. 
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ABSTRACT 

An  experimental-numerical  hybrid  technique  has  been 
developed  to  predict  the  intergranular  contact  load  transfer  in 
granular  media  subjected  to  explosive  loading.  The  granular  media 
was  simulated  by  assemblies  of  circular  discs  in  contact.  The  peak 
contact  load  transfer  coefficients,  i.e.  the  ratio  of  the  maximum 
output  contact  load  to  the  input  contact  load  as  a  function  of  the 
contact  angle,  of  a  given  panicle  were  obtained  through  controlled 
experimental  studies.  These  coefficients  along  with  the  principle 
of  superposition  were  then  used  to  predict  the  peak  contact  loads 
in  several  deterministic  as  well  as  random  assemblies  of  discs. 
Although  the  numerical  scheme  was  rather  simple,  the  predicted 
resuits  compared  favorably  with  the  experimental  data  for  several 
different  assemblies. 


NOMENCLATURE 

C  Peak  contact  load  transfer  coefficient 

D  Diameter  of  a  disc 

f.  Material  fringe  value 

h  Thickness  of  a  disc 

N  Fringe  order 

P.  Input  peak  contact  load 

P,  Output  peak  contact  load 

6,  Contact  angle 

X.  Wave  length 

o;,  o,  Principal  stresses 


INTRODUCTION 

Knowledge  of  wave  propagation  in  granular  materials  is  of 
importance  in  many  branches  of  engineering.  These  include 
powder  metallurgy,  transducer  design,  earthquake  engineering,  soil 
mechanics  etc.  Granular  powders  are  of  great  importance  to  the 
forming  of  many  solid  materials.  These  materials  for  one  reason 
or  another  can  not  be  extruded,  rolled  or  drawn,  molded  and  fired, 
grown  as  crystals  or  cast  from  a  melt,  and  so  are  frequently 


sintered  by  either  static  or  dynamic  compression  of  powders.  Such 
sintering  processes  are  characterized  by  load  transfer  and  flow  of 
panicles  in  a  granular  media.  In  addition,  granular  materials  are 
excellent  shock  wave  attenuators  and  as  such  have  applications  as 
diverse  as  transport  packing  matenals  and  shock  isolation  materials 
for  explosive  chambers.  A  granular  medium  is  composed  of  a 
large  number  of  distinct  particles  as  well  as  some  voids.  The 
particles  can  displace  independently  from  one  another.  The  voids 
may  be  filled  with  gas,  as  in  a  sample  of  dry  sand,  or  they  may 
be  filled  with  liquids,  as  in  a  sample  of  wet  sand.  The  discrete 
character  of  the  granular  media  results  in  a  quite  complex  behavior 
under  conditions  of  dynamic  loading. 

The  study  of  granular  media  from  a  microstructure  point  of 
view  started  many  years  ago  with  research  modelling  granular 
media  as  aggregate  assemblies  of  discs  or  spheres.  Early  work  by 
Iida  [1],  Takahashi  and  Sato  [2],  Hughes  and  Cross  [3],  Gassman 
[4]  and  Brandt  [51  employed  a  normal  granular  contact  force 
concept,  while  Duffy  and  Mindlin  [6]  included  tangential  forces 
as  well.  This  initial  work  investigated  wave  propagation  velocity 
as  a  function  of  confining  pressure,  particle  size  and  aggregate 
geometry.  Important  new  work  has  been  done  by  Oda  [7],  Nemat- 
Nasser  and  Mehrabadi  [8],  et  al.  in  the  development  of  so-called 
fabric  tensor  theories  for  granular  materials.  The  development  of 
distributed  body  theory  by  Goodman  and  Cowin  [9]  has  also 
proved  useful  in  modelling  such  materials.  Computational  work  in 
modeling  large  aggregate  assemblies  of  particles  has  been  carried 
out  by  using  the  distinct  element  method  of  Cundall  and  Strack 
[10],  see  also  Sadd  et  al.  [11]. 

With  regard  to  experimental  techniques,  the  method  of 
photoelasticity  provides  whole  field  data  during  an  experiment  and 
makes  it  possible  to  determine  the  contact  load  between 
neighboring  granular  particles.  Dantu  [12]  and  Wakabayashi  [13], 
Drescher  and  de  Josselin  de  Jong  [14],  and  Durelli  and  Wu  [15] 
have  employed  photoelasticity  to  study  the  granular  media  behavior 
under  static  lolling.  Shukla  et  al.  [16]  [17]  [18]  have  used 
photoelasticity  combined  with  high  speed  photography  to 
investigate  the  wave  propagation  phenomena  due  to  explosive 
loading  of  an  assembly  of  discs.  Experimental  methods  provide 
sufficient  information  to  determine  wave  speeds,  dynamic  contact 
loads,  etc.  in  different  kinds  of  granular  packings  under  various 
loading  conditions. 


DI 


An  experimental-numerical  hybrid  technique  for  analyzing  the 
dynamic  load  transfer  in  a  granular  medium  is  discussed  in  this 
paper.  In  the  experiments,  dynamic  loading  was  achieved  by 
detonating  a  small  cnarge  of  PETS'  in  a  specially  designed  charge 
holder,  which  was  mounted  at  the  top  of  the  experimental  model. 
High  speed  photography  and  dynamic  photoelasticity  were  used 
to  experimentally  determine  the  peak  contact  loads  between  all 
the  panicles.  The  granular  medium  was  simulated  by  assemblies 
of  circular  discs  in  various  packing  geometnes  as  shown  in  Fig.l. 
The  microstructure  of  the  granular  medium  can  be  characterized 
by  branch  vectors  drawn  from  the  mass  centers  of  adjacent  discs 
as  shown  in  Fig.  2.  The  angle  9,  as  shown  in  Fig.  2,  is  drawn 
between  any  two  neighboring  branch  vectors  and  is  denoted  as 
the  contact  angle.  Initially  a  series  of  calibration  experiments  of 
controlled  microstructure  were  conducted  in  which  the  contact 
angles  between  the  granules  were  systematically  varied  The  data 
obtained  from  each  experiment  was  analyzed  to  get  a  load  transfer 
coeffHent,  i.e  a  ratio  of  the  maximum  output  contact  load  to  the 
input  contact  load  for  a  given  disc.  This  load  transfer  coefficient 
was  obtained  for  several  contact  angles  and  number  of  contacts 
per  disc.  With  various  load  transfer  coefficients  determined,  a 
numerical  scheme  was  developed  using  the  experimental  data  and 
based  on  the  principle  of  superposition.  Predictions  were  made  of 
the  intergranular  contact  loadings  for  several  model  assemblies 
with  different  microstructure  as  shown  in  Fig.  1.  Numerical 
predictions  were  compared  with  experimental  data  for  these 
models.  Although  the  numerical  scheme  is  very  simple,  the 
predicted  results  compared  fairly  well  with  the  photoelastic 
experimental  data. 


DETERMINATION  OF  THE 
TRANSFER  COEFFICIENTS 

Granular  media  transmit  mechanical  loadings  primarily  through 
contact  mechanisms  between  adjacent  particles.  This  phenomenon 
is  quite  a  complex  process  and  strongly  depends  on  the  contact 
angle  between  the  particles.  For  simplicity  the  granular  materials 
were  simulated  by  assemblies  of  one  inch  diameter,  1/4  inch  thick 
discs  of  Homalite  100.  The  maximum  number  of  the  contact 
points  for  a  disc  in  this  study  is  6  as  shown  in  the  hexagonal 
packing  in  Fig.  3a.  It  was  shown  in  our  previous  work  (19|  that 
a  mechanical  signal  can  transfer  through  a  contact  point  only  when 
the  contact  angle  is  less  than  90°.  With  respect  to  the  input 
loading  point  of  a  particle,  only  three  contact  points  namely  1,  2 
and  3  as  shown  in  Fig.  3a  will  have  contact  angles  less  than  90°. 
So  if  the  mechanical  signal  enters  a  disc  through  a  contact  point, 
only  three  contact  points  will  be  able  to  transmit  the  response.  In 
this  study,  all  the  particle  interactions  in  the  granular  assembly 
were  represented  by  a  general  four  contact  point  model  as  shown 
in  Fig.  3b.  For  the  case  of  9,  >  90°  and  9,  >  90°  only  two 
contact  points  may  transfer  dynamic  signals,  as  shown  in  Fig.  3c, 
whereas  for  the  case  of  9,  <  90°  and  9,  >  90°,  dynamic  signals 
can  be  transmitted  through  three  contact  points  as  shown  in  Fig. 
3d. 

In  order  to  get  the  peak  contact  load  transfer  coefficients, 
experiments  were  conducted  on  the  three  groups  of  models,  that 
is,  the  two,  three  and  four  contact  point  models  shown  in  Fig  3b, 
3c  and  3d.  In  experiments  of  the  two  contact  point  model,  contact 
angles  9,  and  9,  were  kept  at  90°.  Only  contact  angle  9,  was 
systematically  varied  from  0°  to  90°  In  experiments  of  the  three 
contact  point  model,  contact  angle  9,  was  kept  at  90°,  and  both 
contact  angles  9,  and  9j  were  systematically  varied.  In 
experiments  of  the  four  contact  point  model,  all  three  contact 
angles  9,,  9j  and  9,  were  systematically  varied.  The  experimental 
models  were  placed  in  the  optical  bench  of  a  high  speed  multiple 
spark  gap  camera.  The  camera  was  triggered  at  some  prescribed 


delay  time  after  igniting  the  explosive.  The  high  sptro 
photographic  system  operated  as  a  senes  of  high  intensity 
extremely  short  duranon  pulses  of  light  and  provided 
isochromatic  fringe  images  at  discrete  times  dunng  the  dynamic 
event. 

A  careful  inspection  of  the  photographs  obtained  from  the 
expenments  revealed  that  the  wave  length,  X,  of  the  loading  pulse 
was  much  larger  than  the  disc  diameter  D,  in  fact  X=4D. 
Furthermore,  in  most  cases  the  fringe  patterns  around  the  contact 
points  are  symmetric  on  either  side  of  the  contact  points  and 
similar  to  the  fringes  obtained  under  static  compression.  Both 
these  features  indicate  that  around  the  contact  zone,  quasi-static 
loading  was  present  during  the  wave  propagation  event.  Thus 
Hertz  contact  stress  theory  can  be  used  to  estimate  the  stress  field 
in  the  vicinity  of  the  contact. 

The  isochromatic  fringes  photographed  during  the  expenmeius 
are  lines  of  constant  maximum  shear  stress,  and  are  related  to  the 
stress  field  by  the  stress  optic  law 

o,  -  ct2  =  Nf„  /  h  f  1 ) 

where  o,  and  o2  are  the  principal  stresses,  N  is  the  fringe  ordei 
f„  is  the  material  fringe  value,  and  h  is  the  model  thickness 

In  equation  (1),  o,  and  o2  are  substituted  from  Hertz  contuci  i 
stress  field  equations,  involving  unknowns  of  contact  length  b  and 
friction  factor  (3.  This  equation  was  solved  using  an 
overdeterministic  method  developed  by  Shukla  and  Nigam  [20]  to  1 
accurately  determine  the  contact  length  and  friction  factor  from 
the  full  field  photoelastic  fringe  patterns.  These  obtained  values 
were  substituted  in  the  Hertz  stress  field  equations  and  the  contar-i 
stresses  were  numerically  integrated  along  the  contact  length  u, 
obtain  the  normal  and  tangential  contact  loads. 

Thu  peak  contact  load  transfer  coefficient  C,  at  contact  ponu 
j  is  a  function  of  the  contact  angles  9,,  02  and  9,  and  is  defined 
as  the  ratio  of  the  output  peak  contact  load  P„  to  the  input  peak 
contact  load  P,  ,  i.e. 

C,  =  C,(9,  ,  92 , 9j)  =  P„  /  P,  j=l,  2,  3  (2; 

For  a  two  contact  point  model,  as  shown  in  Fig.  3c,  C,  and  C, 
are  identically  equal  to  zero  since  contact  angles  9,  and  9,  art 
90°,  while  for  three  contact  point  model,  as  shown  in  Fig.  3d,  C, 
is  identically  zero. 

The  transfer  coefficients  C,  thus  obtained  from  the  three  groups 
of  experiments  were  plotted  as  a  function  of  the  contact  angles  9,, 

92  and  9„  as  shown  in  Figs.  4  to  6.  In  Fig.  4,  curve  1  represents 
the  ratio  of  P,j/P,  vs.  contact  angle  02  when  both  the  contact  angles 
9,  and  9,  are  equal  to  90°.  Since  both  P„,  and  P„,  are  zero,  it 
actually  represents  the  transfer  coefficient  Q  of  the  two  contact 
point  model.  The  remaining  curves  in  Fig.  4  represent  the  transfer 
coefficients  Q  vs.  contact  angle  92  of  the  three  contact  point 
model  when  contact  angle  9,  is  90°  and  9,  is  also  a  constant  but 
less  than  90°  (it  is  equal  to  30°,  45°,  60°  and  75°  respectively). 
Fig.5  shows  the  transfer  coef  ents  C,  vs.  contact  angle  9,  for  the 
four  contact  point  model  wnen  the  contact  angles  92=0°  with 
various  values  of  9,,  while  Fig.  6  shows  the  transfer  coefficient  C 
vs.  contact  angle  9,  for  the  same  four  contact  point  model.  The 
transfer  coefficient  C,  for  the  four  contact  point  model  can  be 
obtained  easily  using  Figs.  4  and  6  and  the  property  that 
C,(0„02,9OO)  =  QM.^O0)  and  C,(0„O°,0,)  =  C,(9„O°.0,). 


NUMERICAL  METHOD 

The  experimentally  determined  peak  contact  load  transfer 
coefficients  were  used  to  construct  a  numerical  scheme  capable  of 
predicting  dynamic  load  transfer  in  granular  aggregate  assemblies 
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loading  amplitude.  Final!'.  simple  superposition  of  loading,  as 
mown  in  Fig.  ",  is  used  a. hen  more  than  one  input  contact  load 
occurs  on  a  given  disc. 

Figs.  4  to  6  combined  v*.  uh  a  Lagrangian  interpolation  method 
nave  been  used  to  obtain  all  the  necessary  transfer  c  letticients  to 
predict  the  peak  contact  loads  in  granular  media.  As  an  example 
consider  the  transfer  coefficients  for  a  four  contact  point  model 
: 0  =60' .  9;=0;  and  0,=b5:)  as  shown  in  Fig.  3b.  Since  0  and 
0.  are  equal  to  the  two  contact  angles  of  curve  4  in  Fig.  5.  the 
transfer  coefficient  C,  can  be  obtained  directly  and  is  equal  to 

O. 27 ? .  In  contrast,  the  aansfer  coefficient  C:  can  not  be  obtained 

directly  by  using  Fig.  5.  However  it  can  be  obtained  by  using  the 
curves  in  Fig.  5  combined  with  the  LagTangian  interpolation 
method  and  the  svmmetncal  property  C,(60°.0o.65:’)  = 

C,i65',0°,602.  Using’  the  values  C, (60° ,0260°)  =  028. 
C, (70°  ,0260° i  =  0294'.  Q|8020260°)  =  0.32  and  our  Lagrangian 
interpolation  method,  we  get  C  (6020265°)  =  C,(65'.0260°)  = 
0286  The  transfer  coefficients  Q  can  be  obtained  directly  from 
curve  3  in  Fig.  6.  C.  =  C(6020265°)  =  0.775. 

After  the  three  transfer  coefficients  have  been  obtained,  the 
relevant  output  peak  contact  loads  for  the  four  contact  point  model 
can  be  computed  easily.  According  to  the  definition  of  the  peak 
contact  load  transfer  coefficient,  the  three  output  peak  contact  loads 

P,  ,  P^.  and  Pj,  at  contact  points  1,  2  and  3  are  calculated  as 
follows. 


Po.  =  CjP. 
P„  =  Cd> 


RESULTS  AND  DISCUSSION 

The  experimental- numerical  hybrid  method  was  used  to 
predict  peak  contact  loads  at  every  contact  point  in  various  models 
of  eranular  media.  Four  difterent  microstructural  packings 
illustrated  in  Fig.  1  were  used  in  this  study.  Experimental  fringe 
patterns  obtained  for  each  of  the  microstmctures  are  shown  in  Figs. 

8  to  11.  The  comparison  of  the  numerical  and  experimental 
results  are  shown  in  Figs.  12  to  15. 

Fig.  8  shows  a  sequence  of  eight  photographs  obtained  as 
the  wave  travels  in  a  single  chain  granular  medium  for  the 
geometry  of  Fig.  la.  In  this  geometry  each  panicle  has  two 
contact  points,  hence  only  one  transfer  coefficient  is  needed  to 
model  this  geometry.  This  transfer  coefficient  is  obtained  from 
Fig.  4.  Q  =  Q(9020290°)  =  0.97.  The  peak  contact  loads  at 
each  contact  pomt  can  be  determined  as  follows, 

P.  =  QP 

P,  =  QP.  =  CUP  ,J) 


P„  =  CQP, 

The  comparison  of  numerical  and  experimental  peak  contact 
loads  in  the  single  chain  is  shown  in  Fig.  12.  The  average  peak 
contact  load  error  for  this  model  is  computed  to  be  5%.  The 
results  are  in  very  good  agreement  because  most  of  the 
assumptions  made  in  our  numerical  model  are  satisfied. 

Fig.  9  shows  a  sequence  of  four  photographs  obtained  as  the 


*v av e  trave.v  ,n  a  re\.:ron.ii  cloaca  racking  dlCP:  granu.ar 
medium  tor  tne  geomeirv  of  Fie  lb.  In  this  geometry  each 
panicle  has  mx  contact  point-,  The  major  load  in  this  assembly 
i v  transferred  h\  two  pr.marv  chains.  1  and  2.  shown  in  Fig  lb 
Experimental  data  showed  that  the  tangential  contact  loads  were 
very  close  to  zero  for  the  primary  chains,  hence  they  behaved 
similar  to  the  single  chain  assemblies.  In  the  HCP  model  there  are 
only  two  independent  contact  angles,  0=0,  =  60'  and  0;  =  O', 
with  respect  to  the  input  load.  Thus  only  two  transfer  coefficients 
were  needed  for  the  HCP  model.  From  Figs.  5  and  6,  we  obtain 
the  two  transfer  coefficients,  C;(6020260°)  =  C,(6020260°)=  0.28. 
Q(60°,0°,60cl  =  0.75.  These  coefficients  were  used  to  determine 
the  peak  contact  loads  throughout  the  assembly.  The  comparison 
of  numerical  and  experimental  results  of  the  peak  contact  loads  for 
the  HCP  model  is  shown  in  Fig.  13.  The  average  error  of  the 
peak  contact  loads  for  this  model  is  13%,  with  the  average  error 
along  the  center  line  of  the  primary  chain  being  only  6%. 

Fig.  10  shows  a  sequence  of  four  photographs  obtained  as 
the  wave  travels  in  the  geometry  of  Fig.  lc,  which  will  be  referred 
to  as  a  half  hexagonal  closed  packing  (HHCP)  granular  medium. 
In  this  geometry  a  panicle  has  either  four  or  five  contact  point:-. 
It  is  observed  that  most  of  the  energy  was  transferred  through  a 
vertical  column  consisting  of  the  FIHCP  cells  under  the  explosive 
and  several  horizontal  chains  as  shown  in  Figs.  10  and  14.  The 
peak  contact  loads  were  obtained  by  the  same  method  discussed 
previously.  The  comparison  of  numerical  and  experimental  peak 
contact  loads  is  shown  in  Fig.  14.  It  was  found  that  the  average 
error  for  this  model  is  12%.  However  the  average  error  along  the 
center  line  of  the  horizontal  chains  is  only  7.5%. 

Fig.  1 1  shows  a  sequence  of  four  photographs  as  the  wave 
travels  in  a  random  packing  granular  medium  for  the  geometry  of 
Fig.  Id.  In  this  geometry  particles  have  contact  points  ranging 
from  two  to  six.  The  fringes  in  Fig.  1 1  reveal  a  complex  nature 
of  load  transfer  phenomenon.  The  energy  transfer  showed  no 
preferential  direction  in  this  model.  In  the  former  three  models, 
the  tangential  contact  loads,  especially  along  the  main  path  of  the 
energy  transfer,  were  quite  small.  However  in  this  random  model, 
at  the  contacts  near  the  explosive  point,  the  fringes  appeared 
unsymmetneal  with  respect  to  the  contact  points.  So  it  appears 
that  sizable  tangential  loadings  existed  in  this  case.  Away  from 
the  explosive  point,  the  fringes  showed  the  tendency  to  become 
symmetric  with  respect  to  the  contact  points.  Again  from  Figs.  4 
to  6  all  the  transfer  coefficients  necessary  for  the  determination  of 
the  peak  contact  loads  in  the  random  packing  granular  medium 
were  obtained.  The  comparison  of  numerical  and  experimental 
peak  contact  loads  is  shown  in  Fig.  15.  The  tangential  contact 
loads  tend  to  increase  the  average  error  in  the  peak  contact  loads. 
This  error  was  computed  to  be  19  0%  for  this  model. 


SUMMARY 

A  hybrid  experimental-numerical  technique  has  been 
developed  to  predict  dynamic  contact  loadings  due  to  explosive 
loading  in  different  assemblies  of  circular  discs.  For  a  known 
geometrical  arrangement  of  the  discs  the  technique  can  predict 
contact  loads  at  any  point  in  the  assembly  for  a  given  input 
loading.  The  method  utilizes  experimentally  generated  load 

transfer  coefficients  along  with  simple  linear  superposition  in 
space  The  results  from  this  scheme  are  compared  with  those 
obtained  experimentally  using  the  method  of  dynamic 
photoelasticity.  In  general,  the  results  are  in  good  agreement  for 
regular  packings  of  the  discs.  However  for  a  random  packing  the 
agreement  is  marginal,  and  this  is  primarily  due  to  the  fact  that  the 
numerical  scheme  currently  does  not  take  into  account  tangential 
contact  loads  which  were  quite  large  in  random  arrangement. 
Further,  the  superposition  method  does  not  account  for  any  angular 


D3 


dependence  of  wave  length  of  the  loading  pulse.  It  was  shown  in 
our  previous  work  [191  that  the  duration  of  contact  loads  is 
dependent  on  the  contact  angles.  Thus  to  obtain  better  predictions, 
superpositions  must  be  used  both  in  space  as  well  as  time. 
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ABSTRACT 

Dynamic  load  transfer  in  granular  material  occurs  essentially  through  contact  mechanisms  between  each 
grain.  This  phenomenon  is  quite  a  complex  process  and  depends  inherently  on  the  microstructural  packing 
arrangement  of  the  media.  An  experimental  investigation  is  conducted  to  study  wave  propagation  in 
systematic  aggregates  of  disc  assemblies.  Attention  is  focused  on  load  transfer  path,  wave  velocities,  wave 
attenuation  and  the  stresses  which  are  generated  at  the  contacts  due  to  the  passage  of  stress  waves.  The 
optical  technique  of  photoelasticity  along  with  high  speed  photography  is  utilized  to  get  the  whole  field 
stress  distribution  in  the  granular  media  which  is  subjected  to  explosive  loading.  The  granular  media  is 
simulated  by  circular  discs  made  of  photoelastic  materials.  Several  systematic  arrangements  of  these  discs  are 
used  in  different  experiments.  The  results  show  that  the  load  transfer  phenomenon  depends  strongly  on  the 
rmcrostructure  of  the  media.  Experiments  are  performed  with  initially  unstressed  media  as  well  as  media 
with  prestress. 

1.  INTRODUCTION 

It  is  commonly  observed  that  the  microstructural  packing  or  the  fabric  of  the  granular  media  strongly 
affects  its  mechanical  behavior.  Considerable  effort  has  been  spent  both  experimentally  and  theoretically  to 
look  at  the  influence  of  microstructure  on  c  .formation  and  strength  of  granular  media.  Oda  [1,2]  was  among 
the  first  to  observe  the  microstructure  or  fabric  of  sand  in  unloaded  and  deformed  specimens.  These 
specimens  were  tested  in  the  triaxial  and  direct  shear  apparatus  by  using  thin  sections  and  a  polarized 
microscope.  Borowicka  and  Arthur  [3, 4]  observed  and  measured  the  microstructural  change  of  sand  loading 
by  means  of  either  microscopy  or  radiography.  Konishi  [5]  conducted  biaxial  compression  and  simple  shear 
tests  on  two  dimensional  granular  media,  fabricated  from  epoxy  resin  cylinders,  and  examined  the 
microstructure  change  and  the  state  of  stress  in  the  granules  using  photoelasticity.  All  this  experimental 
work  and  that  of  several  others  [6, 7]  suggest  that  the  concept  of  the  fabric  of  the  granular  media  can  be  very 
useful  in  understanding  their  mechanical  properties.  Based  on  these  experimental  findings,  several  theoretical 
[8,  9]  models  have  developed  constitutive  equations  for  the  granular  media  taking  into  account  the 
microstructural  characteristics.  However,  most  of  this  work  has  been  done  for  static  loading  conditions. 
Recently,  Shukla  and  Damania  [10]  have  conducted  dynamic  experiments  using  photoelasticity  to  study  wave 
propagation  in  granular  media.  Their  experiments  show  that  the  dynamic  load  transfer  depends  on  the  angle 
made  by  the  normals  at  the  contact  points  of  two  adjacent  granules. 

The  purpose  of  this  study  is  to  evaluate  the  effect  of  different  microstructure  or  fabric  of  granular  media 
on  the  wave  propagation  phenomenon.  A  total  of  five  different  regular  microstructural  arrangements  as 
shown  in  Fig.  1  are  considered  for  the  experiments.  The  choice  of  these  microstructural  arrangements  is 
based  on  the  work  of  Shahinpoor  [11]  who  has  demonstrated  that  the  randomly  oriented  granules  follow 
certain  distinct  packing  geometries  reffered  to  as  Vomoi  cells.  High  speed  photography  was  used  to  record 
dynamic  isochromatic  fringes  in  the  birefringen  granules  as  a  function  of  time  due  to  the  passage  of  stress 
waves.  Dynamic  contact  stresses  at  each  conta  ;  lint  and  average  wave  velocities  were  calculated  from  the 
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expei  imental  data.  The  effect  of  microstmcture  on  these  quantities  was  evaluated.  An  experiment  was  also 
conducted  to  look  at  the  effect  of  prestress  on  the  wave  propagation  phenomenon  in  the  body  centered  cubic 
geometry. 

2.  EXPERIMENTAL  PROCEDURE 


The  geometries  of  the  models  used  in  this  study  are  shown  in  Fig.  1.  Five  different  microstructural 
packings  or  cells  based  on  the  work  of  Shahinpoor  [1 1]  were  investigated.  The  granular  media  was  simulated 
with  one  inch  diameter  1/4  inch  thick  Homalite  100  discs.  Homalite  100  is  a  brittle  polyester  material 
which  becomes  temporarily  birefringent  under  the  application  of  load.  Dynamic  loading  was  achieved  by 
detonating  a  small  charge  of  PETN  in  a  specially  designed  charge  holder,  which  was  mounted  either  at  the 
top  or  at  the  center  of  the  model  assemblies  as  also  shown  in  Fig.  1.  The  experimental  models  were  placed 
in  the  optical  bench  of  a  high  speed  multiple  spark  gap  camera.  The  camera  was  triggered  at  some  prescribed 
delay  time  after  igniting  the  explosive.  This  high  speed  photographic  system  operated  as  a  series  of  high 
intensity,  extremely  short  duration  pulses  of  light  and  prrvided  20  photoelastic  images  at  discrete  times 
during  the  dynamic  event  A  typical  sequence  of  four  images  from  each  experiment  are  shown  in  Figs.  2-7. 

3.  ANALYSIS  PROCEDURE 


A  careful  inspection  of  the  photographs  obtained  from  the  experiments  revealed  that  the  wave  length  \ 
of  the  loading  pulse  was  much  larger  than  the  disc  diameter  D(X.  =  4D).  Furthermore,  the  fringe  pattern 
around  the  contact  points  were  symmetric  on  either  side  of  the  contact  points  and  were  similar  to  the  fringes 
obtained  under  static  diametral  compression.  Both  these  features  indicated  that  around  the  contact  zone,  quasi 
static  loading  was  present  during  the  wave  propagation  event  Thus  Hertz  equations  were  used  to  obtain  the 
contact  stresses,  strains  and  loads. 

From  the  Hertz  contact  stress  theory,  the  stress  field  equations  around  the  contact  region  of  two  bodies, 
as  shown  in  Fig.  8,  are  represented  as 
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EXPLOSIVE 


(a) 


(b) 


(c) 


Figure  1 .  Micros  true  turaJ  arrangements  used  in  different  experiments. 
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Figure  7.  Isochromatic  fringes  obtained  in  a  prestressed  model,  geometry  If. 
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Subscripts  1  and  2  refer  to  the  bodies  making  the  contact  as  shown  in  Fig.  8.  Rj,  Rj',  R2  and  R2'  are  the 
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principal  radii  of  curvature  at  the  point  of  contact  and  a  is  the  angle  between  the  corresponding  planes  of  the 
principal  curvatures.  E  is  the  modulus  of  elasticity  and  v  is  Poisson's  ratio. 

The  stress  field  equations  are  coupled  with  the  stress  optic  law 


Nf 


where  aj,  aj  are  the  principat  stresses,  N  is  the  fringe  order,  fa  is  the  material  fringe  value,  and  h  is  the 
model  thickness. 

In  order  to  accurately  determine  the  contact  length  and  friction  factor  from  the  full  field  photoelastic  data, 
Eq.  (4)  is  solved  using  an  overdeterministic  method  developed  by  Shukla  and  Nigam  [12].  These  values  were 
substituted  in  the  Hertz  stress  field  Eqs.  (1  -  3)  and  the  contact  stresses  were  numerically  integrated  along  the 
contact  length  to  obtain  the  contact  loads. 

4.  RESULTS  AND  DISCUSSION 

The  dynamic  isochromatic  photographs  obtained  from  all  the  experiments  showed  the  full  field 
development  of  the  load  transfer  process  in  the  bulk  of  the  granular  media.  These  photographs  clearly  reveal 
the  effect  of  microstructure  on  the  load  transfer  process  and  provided  the  data  necessary  to  obtain  the  wave 
velocity  in  every  direction  and  also  the  normal  and  tangential  load  at  every  contact  point. 

Fig.  2  shows  a  sequence  of  four  photographs  obtained  during  wave  propagation  in  a  single  chain  of  discs 
(geometry  la).  The  wavelength  of  the  pulse  is  seen  to  be  around  four  discs  diameters.  Also  the  fringes  are 
normal  to  the  contact  and  there  are  no  fringes  due  to  the  side  supports  of  the  loading  frame.  This  clearly 
illustrates  that  there  is  no  energy  transfer  occurring  at  the  side  supports  and  that  all  the  energy  is  channelled 
down  the  chain.  From  the  photographs  the  wave  front  position  was  plotted  as  a  function  of  time.  The 
gradient  of  this  graph  gave  the  instantaneous  wave  velocity.  It  was  observed  that  the  wave  velocity  drops 
from  1240  m/s  to  1000  m/s  in  the  initial  four  inches  of  travel.  After  this  the  wave  travels  with  a  constant 
velocity  of  1000  m/s.  The  reason  for  the  initial  high  velocity  may  be  that  the  first  few  discs  are  compressed 
stronger  against  each  other  and  the  ratio  of  the  incremental  strains  to  incremental  stresses  might  be  very 
large.  This  will  lead  to  larger  effective  modulus  and  thus  higher  velocity. 

The  contact  loads  obtained  from  photoelastic  fringes  were  plotted  for  different  contacts  as  a  function  of 
time  as  shown  in  Fig.  9.  At  a  given  contact  the  load  increased  from  zero  to  a  peak  value  and  then  gradually 
decreased  to  zero.  A  typical  duration  of  the  pulse  at  the  contact  was  1 10  ps.  The  wavelength  of  the  pulse 
was  calculated  by  taking  the  product  of  average  wave  speed  and  the  pulse  duration  at  the  contact.  The 
wavelength  thus  obtained  was  120  mm  (i.  e.,  around  four  disc  diameters).  This  value  compared  well  with  the 
visual  observation  from  the  photographs.  Due  to  internal  losses  within  the  granule,  energy  spent  in  closing 
the  contact  and  some  frictional  and  reflection  effects,  the  peak  contact  loads  dropped  as  the  distance  from  the 
explosive  loading  increased.  Using  the  values  of  peak  loads,  it  was  found  that  there  was  a  20%  decrease  in 
load  as  the  wave  travelled  5  disc  diameters  starting  from  the  second  disc.  This  is  much  higher  when 
compared  with  the  drop  in  peak  load  for  a  uniform  bar  which  is  2%  for  the  same  distance  of  travel. 

Fig.  3  shows  a  sequence  of  four  photographs  obtained  due  to  wave  propagation  in  geometry  lb.  This 
geometry  has  a  coordination  number  of  three,  i.e.,  each  grain  is  in  contact  with  three  other  grains.  The 
wavelength  of  the  loading  pulse  is  about  5  disc  diameters.  Again,  no  fringes  appear  at  the  contacts  with  the 
supporting  loading  frame  indicating  that  all  the  energy  is  channelled  along  the  cell  structure.  The  average 
wave  velocity  in  the  vertical  direction  was  found  to  be  8G0  m/s.  This  is  about  20  percent  lower  than  the 
average  velocity  in  the  single  chain  disc  assembly.  This  decrease  in  the  average  velocity  is  largely  due  to  the 
fact  that  the  wave  velocity  is  largest  in  the  direction  normal  to  the  contact  points.  When  contact  points 
between  two  adjacent  granules  deviate  from  180°  this  velocity  drops.  Thus  the  wave  velocity  in  this 
assembly  is  higher  when  the  transfer  occurs  from  a  granule  which  is  vertically  on  top  of  the  other,  whereas 
the  velocity  drops  when  the  granules  are  at  an  angle  to  each  other. 

Besides  the  wave  velocity,  the  photoelastic  fringes  were  also  analyzed  to  get  the  contact  loads  as  a 
function  of  time  as  shown  in  Fig.  10.  The  peak  contact  load  at  the  contact  point  1  was  found  to  be  BOON. 
As  the  wave  entered  the  following  two  contacting  granules  the  peak  load  dropped  by  40%  at  contact  points  2 
and  3.  The  loading  wave  then  enters  the  next  granule  through  contact  points  4  and  5.  Since  these  points  are 
not  normal  to  points  2  and  3  the  peak  loads  further  drop  by  35  percent.  The  contact  load  at  point  4  builds  up 
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as  energy  which  is  channelled  in  through  two  contact  points  has  only  one  exit  point.  The  peak  load  at  point 
6  wa$  found  to  be  700N.  Thus  in  going  from  point  1  to  6,  a  vertical  distance  equal  to  2.73  times  the 
diameter  of  the  disc,  there  is  an  attenuation  of  about  46  percent  in  the  peak  load.  This  is  considerably  larger 
than  the  drop  seen  in  the  single  chain  assembly  and  clearly  indicates  the  effect  of  microstructure  on  the  wave 
propagation  phenomenon. 

The  wave  propagation  phenomenon  in  a  Body  Centered  Cubic  (BCC)  arrangement,  geometry  lc,  is 
shown  in  the  four  frames  of  Fig.  4.  The  wave  propagated  only  in  a  single  chain  in  this  arrangement.  No 
energy  was  transferred  to  other  chains  in  the  assembly.  Comparison  of  Fig.  4  and  Fig.  2  shows  that  the 
wave  propagation  phenomenon  is  identical  for  the  single  chain  and  the  BCC  packing,  when  the  explosive 
loading  is  on  the  top  of  a  disc  in  the  assembly. 

The  photoelastic  fringes  for  the  hexagonal  closed-packing  arrangement,  geometry  Id,  are  shown  in 
Fig.  5.  The  wave  propagation  in  this  experiment  does  occur  in  two  dimensions.  The  load  transfer  in  this 
experiment  can  be  categorized  by  two  distinct  chains:  the  primary  chains  (such  as  chains  a  and  b  as  shown  in 
Fig,  Id)  which  emanate  from  the  disc  on  which  the  explosive  loading  takes  place  and  the  secondary  chains 
which  emanate  due  tc  contact  of  other  discs  with  the  discs  in  the  primary  chain. 

The  normal  and  tangential  loads  were  computed  from  the  photoelastic  data  at  each  point  as  a  function  of 
time.  It  was  observed  that  the  tangential  loads  were  very  close  to  zero  for  the  primary  chains.  Thus,  the 
contacts  of  the  primary  chains  were  frictionless,  indicating  that  they  behaved  similar  to  the  single  chain 
assemblies  but  with  inclined  orientation  with  respect  to  the  explosive  loading.  The  contact  loads  versus  time 
for  different  contacts  were  plotted  for  the  primary  chains  as  shown  in  Fig.  1 1 .  The  drop  in  the  peak  loads 
from  one  contact  to  another  was  large  as  compared  to  single  chain  or  BCC  packing  experiments,  since  energy 
transfer  occurs  to  the  secondary  chains  as  the  wave  propagates  in  the  primary  chain.  The  wave  attenuation 
for  5  disc  diameters  of  travel  in  the  primary  chain  was  70%  as  compared  to  20%  for  the  single  chain  and  the 
BCC  packing.  Also,  the  average  velocity  in  the  vertical  direction  was  950  m/s  which  was  about  the  same  as 
in  a  single  chain  experiment. 

The  peak  normal  contact  load  at  different  contact  points  and  the  load  transfer  paths  are  shown  in  Fig.  12. 
The  peak  loads  in  the  triangular  region  enclosed  by  the  main  chains  are  considerably  higher  than  the  contact 
loads  outside  this  area.  This  indicates  that  the  region  of  intense  loading  due  to  explosion  is  in  the  envelope 
enclosed  by  the  two  primary  chains. 

From  the  BCC  and  HCP  experiments,  it  was  observed  that  the  load  transfer  path  and  the  magnitude  of 
load  transfer  from  one  contact  to  another  was  related  to  the  angle  made  by  the  normals  to  the  contacting  discs 
at  the  contact  point.  If  the  input  loading  is  normal  and  the  angle  made  by  the  normals  between  the  two  discs 
as  shown  in  Fig.  13  is  acute  (ay  <  (rt/2)),  no  load  transfer  would  occur,  however,  if  the  angle  made  by  the 
normals  is  obtuse  (ay  >  (jr/2)),  mechanical  signals  would  propagate  through  the  contact.  To  further 
illustrate  this  point  consider  the  load  transfer  into  and  from  disc  A  in  the  HCP  geometry.  The  normal 
contact  loads  at  every  contact  point  made  by  disc  A  with  the  neighbouring  discs  are  shown  in  Fig.  14.  Wave 
energy  entering  at  point  1  shows  maximum  transfer  across  point  2  which  is  directly  ahead  of  point  1. 
Although  contact  points  3  and  4  make  the  same  angles  with  the  main  chain,  the  contact  loads  are  higher  at 
point  4  than  point  3  due  to  the  superposition  of  loads  from  the  other  discs.  No  load  transfer  occurs  across 
points  5  and  6  as  ay  <  (rt/2)  for  these  contact  points.  A  detailed  study  of  the  angular  dependence  of  load 
transfer  in  a  granular  media  has  been  conducted  by  the  authors  and  is  presented  elsewhere  [13]. 

Fig.  6  shows  a  sequence  of  four  photographs  as  the  wave  travels  in  geometry  le.  This  geometry  has 
coordination  numbers  of  both  4  and  5.  The  isochromatic  fringes  show  a  complex  nature  of  load  transfer 
phenomenon.  The  peak  normal  contact  loads  at  different  contact  points  are  shown  in  the  left  part  of  Fig.  15, 
while  the  main  load  transfer  paths  are  shown  in  the  right.  In  the  vertical  direction  the  bulk  of  the  energy  is 
channelled  through  a  column  consisting  of  the  primary  cell  structure.  In  the  horizontal  direction  there  are 
two  chains  of  load  transfer.  The  first  horizontal  chain  making  contact  with  the  primary  cell  carries  most  of 
the  load  and  is  shown  in  Fig.  15  as  the  high  load  chain.  The  first  horizontal  chain  making  contact  with  the 
second  cell  carries  a  lesser  amount  of  load  and  is  shown  in  Fig.  15  as  the  low  load  chain.  It  appears  that  this 
sequence  is  repeated  as  one  goes  down  the  column  of  the  cells.  The  photograph  from  this  experiment  also 
shows  that  the  energy  entered  each  cell  within  the  column  through  four  contact  points,  ij  to  i4  as  shown  in 
Fig.  15  and  exits  mostly  through  six  contact  points  0j  to  0^.  During  wave  propagation  in  this  geometry, 
position  of  the  maximum  vertical  contact  load  alternated  from  one  cell  to  the  other  in  the  column  between 
the  point  along  the  centerline  (point  A)  and  points  on  the  side  of  the  centerline  (points  B  and  C).  When  the 
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maximum  vertical  contact  load  was  at  the  point  along  the  center  line,  high  load  transfer  occurred  through  the 
horizontal  chain  in  the  adjacent  cells.  However,  when  the  maximum  vertical  contact  load  was  on  the  side  of 
the  centerline  a  lesser  load  was  transferred  through  the  horizontal  chains  in  the  adjacent  cells.  The  average 
wave  velocity  along  the  vertical  direction  in  this  geometry  was  850  m/s  whereas  the  horizontal  chains 
showed  a  velocity  of  1000  m/s.  It  should  be  pointed  out  here  that  when  the  input  load  to  a  disc  had  a  large 
tangential  component,  some  energy  transfer  did  occur  across  to  discs  which  were  at  right  angles. 
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Figure  12.  Peak  normal  contact  loads  a  different  contact  points  and  the  load  transfer  paths  in  the  HCP 
geometry. 


Figure  13. 


Contact  angle  made  by  normals  drawn  from  the  center  of  granules  to  the  contact  point. 
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Finally,  an  experiment  was  conducted  to  study  wave  propagation  in  a  BCC  geometry  with  a  prestress  on 
it.  The  reason  for  conducting  this  experiment  was  to  verify  the  results  of  Fig.  4.  It  was  felt  that  maybe 
wave  propagation  in  Fig.  4  was  seen  only  in  the  vertical  direction  due  to  gravitational  loading.  Thus  a 
biaxial  prestress  of  250N  was  applied  to  the  BCC  arrangement  of  discs  using  a  specially  designed  loading 
fixture  as  shown  in  Fig.  16.  The  dynamic  explosive  loading  was  applied  in  the  center.  A  series  of  four 
photographs  obtained  from  this  experiment  are  shown  in  Fig.  7.  They  verifiy  that  wave  propagation  only 
occurs  along  the  two  normal  chains  with  no  energy  transfer  across  contacts  which  are  at  right  angles.  The 
only  difference  in  this  experiment  was  that  the  velocity  was  about  25%  higher  in  comparison  to  the 
experiment  without  any  prestress. 


Figure  16.  Biaxial  loading  frame. 


5.  SUMMARY 

Dynamic  photoelasticity  was  employed  to  study  the  effect  of  microstructure  on  wave  propagation  and 
dynamic  load  transfer  in  a  granular  media.  The  granular  media  was  modeled  as  one  and  two  dimensional 
arrays  of  circular  discs  fabricated  from  photoelastic  material  Homolite-100.  The  experimental  data  were 
analyzed  to  determine  the  wave  velocities,  identify  characteristic  dynamic  load  transfer  paths  and 
quantitatively  calculate  the  dynamic  contact  forces  at  each  contact  point. 
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An  experimental  investigation  was  conducted  to  study  the  dependence  of  packing  geometry  on  the  dynamic  load 
transfer  in  two  dimensional  granular  aggregate  chains.  The  granular  media  was  simulated  by  circular  discs  made  of 
photoelastic  material.  The  experimental  method  utilizes  the  combination  of  high  speed  photography  and  photoelasu- 
city  to  visualize  the  dynamic  load  transfer  phenomenon  in  two  dimensional  model  chains  subjected  to  explosive 
loading  The  photographs  thus  obtained  were  analysed  to  get  the  normal  and  tangential  loads  at  the  contact  points  as 
a  function  of  time.  The  results  indicate  that  the  load  transfer  is  strongly  dependent  on  the  angle  between  the  vectors 
drawn  from  the  mass  centres  of  the  contacting  granules. 


1  INTRODUCTION 

Granular  media  transmit  mechanical  loadings  primarily 
through  contact  mechanisms  between  each  grain.  This 
phenomenon  is  quite  a  complex  process  and  depends 
inherently  on  the  microstructural  packing  arrangements 
of  the  media.  Porosity  alone  is  not  a  sufficient  measure  to 
characterize  such  a  load  transfer  process.  Recently, 
several  theoretical  and  experimental  investigations 
(lX2)t  have  been  conducted  to  relate  microstructure  to 
macroscopic  behaviour.  A  general  finding  of  some  of  this 
work  is  that  local  microstructure  or  fabric  is  significant 
and  that  particular  fabric  vectors  and  tensors  can  be 
used  to  develop  theories  to  predict  the  mechanical 
behaviour  of  such  materials.  In  particular,  branch 
vectors  between  the  mass  centres  of  typical  grains  and 
normal  vectors  in  the  direction  of  the  contact  normals 
have  been  proposed  (see  Fig.  1).  Specifically,  Nemat- 
Nasser  et  al.  (1)  have  suggested  mechanical  constitutive 
relations  based  upon  writing  the  stress  as  a  function  of  a 
second  order  fabric  tensor,  Ft),  where 


with  M  being  the  number  of  n-contacts  per  unit  volume. 
This  study  addresses  this  issue  for  the  case  of  dynamic 
load  transfer  by  investigating  the  effects  of  the  angle 
between  branch  vectors  on  the  wave  propagation 
through  granular  aggregate  assembly  chains. 

The  term  ‘contact  angle'  is  defined  as  the  angle 
between  any  two  branch  vectors  connecting  the  mass 
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centres  of  a  pair  of  granules,  as  shown  in  Fig.  1.  This 
angle  plays  an  important  part  in  determining  the 
dynamic  load  transfer  in  such  a  medium.  Examples  of 
this  fact  have  been  shown  in  the  work  of  Shukla  and 
Damania  (3)  and  are  shown  in  Figs  2  and  3.  These 
dynamic  photoelastic  figures  illustrate  how  waves  from 
an  explosive  charge,  located  at  the  top  of  each  photogra¬ 
ph,  move  through  two  different  granular  packing 
arrangements. 

Studies  of  the  load  transfer  in  granular  media  have 
been  previously  conducted.  Drescher  and  De  Josselin  De 


Fig.  I  Branch  vectors  connecting  the  mass  centres  of  a  pair  of 
granules 
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Frame  3.  I  -  55.5 


Frame  14.  i  =  133  ns 


Fig.  2.  Isochromatic  fringes  due  to 


Jong  (4)  have  simulated  granular  media  by  using 
assemblies  of  circular  discs,  and  then  studied  the  static 
load  transfer  through  the  assembly  by  means  of  photo¬ 
mechanics.  Rossmanith  and  Shukla  (5)  extended  this  idea 
to  the  dynamic  case  through  the  u^e  of  high  speed  pho¬ 
tography.  Additional  dynamic  woric  was  also  carried  out 
by  Shukla  and  Damania  (3),  and  Shukla  and  Rossmanith 
(6).  Most  of  this  previous  work  focussed  on  wave  propa¬ 
gation  phenomenon  in  general  without  going  into  details 
of  relating  specific  microstructure  to  the  associated  wave 
motion. 

The  purpose  of  this  investigation  was  to  study  the 
relation  between  the  dynamic  load  transfer  and  the 
contact  angle  in  single  disc  chains  of  circular  granules. 
This  specific  relation  will  prove  to  be  useful  in  character¬ 
izing  how  local  microstructure  affects  the  wave  propaga- 
tional  phenomena.  Dynamic  photoelasticity  along  with 
high  speed  photography  are  employed  to  collect  time- 
dependent  data  on  the  rapidly  moving  wave  motion. 
Simplified  Hertz  contact  stress  theory  along  with  photo¬ 
mechanics  has  been  used  to  determine  the  load  transfer 
between  pairs  of  granules.  The  wave  motion  of  interest 
here  is  transient  in  nature  being  produced  by  explosive 


Frame  9.  /  =  93  5 


Frame  19.  r  =  164  5  jjs 


in  a  body-centred  cubic  arrangement 


loadings  of  short  duration  yielding  a  primary  wave 
length  of  approximately  four  to  five  grain  diameters. 

2  EXPERIMENTAL  PROCEDURE 

The  experimental  model  used  in  this  study  was  com¬ 
prised  of  disc  chain  assemblies  of  Homalite  100  discs  as 
shown  in  Fig.  4.  Homalite  100  is  a  birefringent  brittle 
polyester  material  whose  mechanical  and  optical  proper¬ 
ties  are  well  characterized  (£  =  4.8  GPa,  v  =  0.35  and 
f„  =  21.9  kN/m).  The  discs  were  1  in  in  diameter  and  1  4 
in  thick.  The  angles  6X  and  92  in  the  assembly  were 
changed  in  a  series  of  experiments.  During  the  expen- 
.ments  the  discs  were  dynamically  loaded  by  firing  a  small 
charge  of  explosive  PETN  in  a  specially  designed  charge 
holder.  The  wave  propagation  phenomenon  due  to 
explosive  loading  in  the  granular  media  was  studied 
using  the  technique  of  dynamic  photoelasticitv  and  high 
speed  photography. 

The  models  were  placed  in  the  optical  bench  of  a  high 
speed  multiple  spark  gap  camera  T  '  e  camera  was  trig¬ 
gered  at  some  prescribed  delay  time  after  igniting  the 
explosive.  This  high  speed  photographic  system  operated 
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Fig.  3.  Isochromatic  fringes  due  to  explosive  loading  in  a  hexagonal  cubic  arrangement 
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as  a  series  of  high  intensity,  extremely  short  duration 
pulses  of  light  and  provided  20  photoelastic  images  at 
discrete  times  during  the  dynamic  event.  A  typical 
sequence  of  three  images  for  two  different  orientations  is 
shown  in  Figs  5  and  6.  These  photographs  of  the  wave 
propagation  process  at  different  stages  of  development 
provided  the  data  necessary  to  obtain  the  contact  loads 
along  different  chains. 


3  ANALYSIS  PROCEDURE 

A  careful  inspection  of  the  photographs  obtained  from 
the  experiments  revealed  that  the  wave  length,  X,  of  the 
loading  pulse  was  much  larger  than  the  disc  diameter,  D. 
The  wave  length  was  determined  by  measuring  the 
length  of  the  photoelastic  fringe  patterns  of  the  loading 
pulse,  and  the  data  indicated  that  x  4 D.  Furthermore, 
the  fringe  patterns  around  the  contact  points  were  sym¬ 
metric  on  either  side  of  the  contact  poir.’s  and  were 
similar  to  the  fringes  obtained  under  static  diametral 
compression  as  shown  in  Fig.  7.  Both  these  features  indi¬ 
cated  that  around  the  contact  zone,  quasi  static  loading 
was  present  during  the  wave  propagation  event.  Thus 


Hertz  equations  were  used  to  obtain  the  contact  stresses, 
strains,  and  loads. 

From  the  Hertz  contact  stress  theory  (7),  the  stress 
field  equations  around  the  contact  region  of  two  discs,  as 
shown  in  Fig.  8,  are  represented  as 
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Fig.  4.  Geometrical  arrangement  used  in  experiments 


where  are 

_ MA  +  B) _ 

(pl  ~  ABsi{2AB  +  2x2  +  2z2  -  2b2) 

_ k(A  -  B) _ 

<t>1  ~  ABj(2AB  +  2x 2  +  2z2  -  2b2) 

A  =  +  *)2  +  z2}< 

B  =  J{(b  -  x)2  +  z2};  A  =  2R  1  —  (3) 

with  E  =  modulus  of  elasticity,  v  =  Poisson’s  ratio,  /3  is  a 
friction  factor,  and  b  is  the  half  contact  length. 

The  stress  field  equations  are  coupled  with  the  stress 
optic  law 


where  <r a2  are  the  principal  stresses,  N  is  the  fringe 
order,/,  is  the  material  fringe  value,  and  h  is  the  model 
thickness. 

In  order  to  accurately  determine  the  contact  length 
and  friction  factor  fr  m  the  full  field  photoelastic  data, 
equation  (4)  is  solved  using  an  overdeterministic  method 
developed  by  Shukla  and  Nigam  (8).  These  values  were 
substituted  in  the  Hertz  stress  field  equations  (2)  and  (3) 
and  'e  contact  stresses  were  numerically  integrated 
alo  te  contact  length  to  obtain  the  contact  loads. 


4  RESULTS  AND  DISCUSSION 
A  series  of  five  groups  of  experiments  was  conducted 
with  the  geometry  shown  in  Fig.  4.  These  groups  include 


Fig.  5.  Typical  isochromatic  fringes  obtained  in  symmetric  arrange¬ 
ment  (0,  =  02) 


a  parametric  study  of  the  effects  of  the  branch  angles  0, 
and  d2  on  load  transfer  phenomenon.  The  branch  angles 
included  the  values  of  30.  45,  60.  75,  and  90  degrees.  The 
dynamic  isochromatic  photographs  obtained  for  two  of 
these  experiments  are  shown  in  Figs  5  a  d  6.  Using  the 
stress  field  equations  (2)  and  (3)  along  with  (4).  and  the 
photoelastic  fringes  near  the  contact  points,  gives  the 
intensity  of  the  contact  stresses.  The  fringes  are  sym- 


JOL'RNAL  OF  STRAIN  ANALYSIS  VOL  :.l  NO  1  IWX  ,  IMcchE  ISHK 


DYNAMIC  LOAD  TRANSFER  IN  TWO-DIMENSIONAL  GRANULAR  AGGREGATES 


metric  in  both  chains,  indicating  predominantly  normal 
loading.  The  absence  of  any  fringe  patterns  at  contact 
points  along  the  side  supports  confirms  that  all  energy  is 
channelled  along  the  disc  chains  and  no  loss  of  energy 
occurs  to  the  side  supports.  The  normal  contact  loads 
obtained  from  these  photographs  for  a  specific  experi¬ 
ment  are  shown  in  Fig.  9  for  case  of  9X  =  75  degrees  and 
0 2  =  30  degrees.  The  three  contact  loads  at  each  of  the 
three  contact  points  are  labelled  P,,  P2,  and  P3  as 
shown.  The  load  at  each  contact  point  increases  as,the 
wave  interacts  with  it,  builds  up  to  a  peak  value,  and 
then  monotonically  decays  to  zero.  It  can  clearly  be  seen 
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that  as  the  branch  angle  increases  the  dynamic  peak  load 
drops.  However,  the  wavelength  of  the  loading  pulse 
increases  with  the  branch  angles.  Figure  10  shows 
another  typical  case  with  9X  =  90  degrees  and  92  =  30 
degrees.  For  this  case  the  load  P3  in  the  90  degree  branch 
was  found  to  be  zero. 

Since  the  significant  point  in  our  study  is  not  the 
actual  loadings  but  rather  the  transfer  characteristics,  the 
ratio  of  transmitted  to  incident  loads  is  to  be  used.  This 
ratio  has  symmetry  with  respect  to  the  branch  angles  9X  , 
and  f?2,  i.e.,  P2/P i  as  a  function  of  9X  and  92  is  the  same 
as  Pj/Pt  written  as  a  function  off?,  and  9l.  Consequently 
only  one  ratio,  P2/Px,  will  be  considered  in  detail. 
Figure  1 1  shows  this  load  transfer  ratio  as  a  function  of 
contact  angles  9X  and  92 .  This  figure  illustrates  the  rapid 
attenuation  of  the  load  transfer  with  the  branch  angle, 
and  it  also  shows  the  inter-relationships  of  (?,  and  (?,  on 
the  load  transfer.  For  example,  a  higher  value  of  f?,  pro¬ 
duces  less  load  transfer  attenuation  for  a  given  branch 
angle,  92 . 

The  effect  of  branch  angles  on  the  signal  wave  length  is 
shown  in  Fig.  12.  Here  the  duration  of  contact  is  plotted 
against  the  angle  92 .  It  is  observed  that  a  significant 
increase  in  this  duration  time  occurs  as  the  branch  angle, 

92 ,  is  increased.  This  means  that  the  granular  assembly 
will  act  as  a  wave  guide  which  will  increase  the  wave 
length  of  the  transmitted  transient  signal. 


5  CONCLUSIONS 

The  experimental  study  conducted  in  this  paper  demon¬ 
strates  the  angular  dependence  of  dynamic  load  transfer 
in  two  dimensional  granular  media.  The  results  indicate 
the  following. 

(1)  Rapid  attenuation  of  load  transfer  occurs  as  the 
branch  angle  increases  from  0  to  90  degrees. 

(2)  The  attenuation  of  load  transfer  also  depends  on  the 
inter-relationship  between  the  two  branch  angles. 

(3)  The  wavelength  of  the  loading  pulse  increases  with 
the  branch  angle.  Here  again  the  inter-relationships 
between  the  branch  angles  has  an  influence  on  the 
duration  of  the  contact  load.  The  wavelength  almost 
doubles  as  the  contact  angle  increases  from  15  to  75 
degrees. 

It  should  be  pointed  out  that  the  current  work  rep¬ 
resents  only  a  first  step  in  understanding  the  wave  propa¬ 
gation  process  in  a  complex  aggregate  assembly.  For 
example,  in  considering  the  wave  motion  in  Fig.  3,  many 
waves  take  very  complex  paths  during  their  propagation 
histories.  The  simple  two  branch  geometries  considered 
in  the  current  study  can  then  be  interpreted  as  only  a 
beginning  in  the  understanding  of  how  the  local  micro- 
structure  effects  the  transmission  of  waves  in  granular 
materia'-..  Finally,  theoretical-numerical  work  is  also 
underway  to  calculate  this  dynamic  load  transfer  pheno¬ 
mena,  and  this  will  be  reported  in  the  near  future. 
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ABSTRACT 

An  experimental  investigation  is  conducted  using  dynamic  photoelasticity  and  high 
speed  photography  to  study  the  wave  propagation  due  to  blast  loading  in  porous  media  as  a 
function  of  fluid  saturation.  The  porous  media  has  been  modelled  as  continuous  solid 
containing  particular  arrays  of  holes  or  voids.  The  study  has  focused  mainly  on  the  effect 
of  the  porous  structure  on  transient  pulse  propagation  as  well  as  the  effect  of  the  moisture 
in  the  pores  on  wave  propagation.  A  series  of  experiments  have  been  conducted  using  a 
sheet  of  Homalite  100  with  different  geometry  of  the  periodic  array  of  holes.  A  small 
amount  of  explosive  is  used  to  generate  the  au^ss  wave.  Dynamic  photoelastic 
photographs  are  taken  with  the  high  speed  camera  as  the  wave  propagates  across  the  holes. 
These  data  are  analyzed  to  obtain  the  wave  velocity  as  well  as  the  stress  wave  attenuation 
in  the  porous  media. 

INTRODUCTION 

The  problem  of  interaction  of .  elastic  waves  with  discontinuities  or  boundaries  of 
complex  shapes  arises  in  situations  where  waves  propagate  through  a  medium  having 
cavities,  inclusions  or  cracks.  Due  to  material  inelasticity  and  inhomogeneity,  the  wave 
propagation  in  a  discontinuous  medium  is  much  more  involved  than  homogenous  elastic 
wave  propagation  and  it  shows  directional  as  well  as  frequency  dependence.  This 
phenomenon  becomes  significant  for  step  loading  pulses  where  the  wavelength  are  of  the 
order  of  the  size  of  the  discontinuities.  Such  problems  denoted  as  scattering  and 
diffraction  problems  have  long  standing  interest  in  acoustics  and  electromagnetic  wave 
theory. 
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Composite  materials  such  as  concrete,  ceramics,  etc.  are  characterized  by  the  number 
of  pores,  voids  or  fluid  filled  cavities.  The  influence  of  these  pores  on  the  deformation 
and  the  failure  of  these  materials  has  not  been  interpreted  uniquely.  Since  areas  of  stress 
concentration  may  arise  in  the  vicinity  of  these  pores,  it  is  believed  that  the  pores  may 
play  a  role  in  influencing  the  crack  and  wave  propagation  in  these  materials.  Moreover, 
porous  materials  are  used  extensively  for  shock  isolation  as  they  are  capable  of  absorbing 
large  quantities  of  energy  during  impact  loading.  Hence  the  behavior  of  these  materials 
under  impulsive  loading  has  been  of  substantial  interest  to  engineers. 

Wave  propagation  in  a  discontinuous  media  has  also  been  of  interest  to  the  soil  and 
rock  mechanics  community.  The  propagation  of  elastic  waves  in  the  earth’s  crust  is  most 
intimately  related  to  the  properties  of  soil  and  rock.  The  elastic  properties  of  these 
substances  are  greatly  affected  by  the  amount  of  water  contained  in  them,  packing  density, 
porosity,  the  size  of  the  particles  that  form  the  substances,  and  the  binding  material  which 
they  contain.  Current  interest  in  geomechanics  is  focussed  on  the  transient  phenomena 
occurring  in  earthquakes,  wave  loading  and  consolidation.  Moreover  the  increasing  needs 
for  urban  and  resource  development  demand  faster,  safer  and  more  efficient  procedures  for 
underground  excavations  of  rock.  Most  methods  of  rapid  excavation  in  hard  rock  use  some 
form  of  dynamic  loading,  such  as  explosive  or  water  jet.  This  type  of  loading  produces 
stress  waves  which  induce  crack  initiation  and  propagation. 

The  initial  attempts  to  study  rock  media  and  soil  structure  as  arrays  of  elastic 
particles  (eg.  spheres  and  discs)  were  made  by  Iida  (1,2),  Takashashi  and  Sato  (3), 
Gasmann  (4)  and  Brandt  (5).  They  investigated  the  propagation  velocity  as  a  function  of 
confining  pressure,  particle  size  and  aggregate  geometry.  The  effect  of  water  content  in 
the  pores  on  elastic  velocity  has  been  studied  by  several  investigators.  Oliphant  (6)  and 
Owen  (7)  found  that  slight  additions  of  w^.ter  caused  a  sharp  drop  in  velocity  with  a  slow 
decrease  as  the  saturation  approached  100  percent.  Hughes  and  Jones  (8)  measured  the 
dilation  wave  velocity  of  samples  of  very  low  porosity,  less  than  1%.  Using  the  same 
apparatus  and  methods,  Hughes  and  Cross  (9)  measured  the  velocity  in  Soienhofen 
Limestone  (porosity  4%)  and  Caplen  Dome  Sandstone  (porosity  5%)  for  dry  and  saturated 
samples. 
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A  considerable  amount  of  research  is  also  under  way  in  determining  the  internal 
structure  of  porous  and  granular  media  by  various  sounding  techniques.  For  example, 
acoustic  emission  methods  have  been  presented  by  Hardy  (10),  while  Aliison  and  Lama 
(11)  have  discussed  a  low  amplitude  vibration  technique  to  predict  rock  structure. 

Current  research  in  wave  propagation  in  such  media  has  involved,  for  example, 
statistical  theories  by  Varadan  et  al  (12)  or  mixture  theories  by  Junger  (13).  Analytical 
approaches  are  often  limited  as  they  cannot  fully  account  for  the  material  inhomogeneity, 
isotropy  and  defects.  Most  of  the  previous  work  focused  on  wave  propagation  phenomena 
in  general  without  going  into  the  details  of  relating  specific  microstructure  to  the  associated 
wave  motion. 

This  paper  reports  on  an  experimental  study  of  wave  propagation  due  to  explosive 
loading  in  a  porous  medium.  The  porous  medium  was  modelled  as  an  array  of  holes 
machined  in  a  continuous  meet  of  a  brittle  polyester  material  Homalite  100.  The  study 
looked  at  the  wave  propagation  phenomenon  from  a  microscopic  point  of  view  by  going 
into  the  details  of  the  geometric  nature  of  the  porous  structure.  The  geometry  of  the  pores 
was  changed  by  varying  the  size  of  the  holes,  changing  the  pitch  or  the  spacing  between 
the  holes  and  changing  the  geometric  arrangement  of  the  holes.  In  all  the  experiments  the 
stress  waves  were  produced  by  explosive  loadings  of  short  duration.  Dynamic 
photoelasticity  and  h;gh  speed  photography  were  employed  to  collect  the  time  dependent 
data  as  the  wave  front  moved  rapidly  through  the  porous  structure. 

EXPERIMENT  AL  PROCEDURE 

To  investigate  the  behavior  of  stress-waves  in  the  periodically  flawed  half-plane,  also 
referred  to  as  the  porous  medium,  a  series  of  dynamic  photoelastic  experiments  was 
conducted.  The  photoelastic  models  were  fabricated  from  a  large  sheet  of  Homalite- 100 
having  dimensions  of  254mm  x  305mm  x  6.4mm.  Homalite-100  is  a  birefringent  brittle 
polyester  material  whose  mechanical  and  optical  properties  arc  well  characterized.  The 
stress-sr  ain  behavior  of  this  material  is  similar  to  that  of  rock.  An  array  of  holes  having  a 
pitch  p  and  diameter  d  was  machined  to  simulate  the  porous  media. 

Figure  1  shows  the  geometry  of  the  porous  model.  Dynamic  loading  to  pm'dnce  the 
incident  wave  was  achieved  by  detonating  80mg  of  lead  azide  directly  on  top  of  the  model. 


G3 


Two  different  arrangements  of  the  holes  were  used  as  shown  in  Fig. 2  The  arrangement  of 
holes  in  the  first  and  second  type  of  configurations  are  referred  to  as  geometry  A  and  B 
respectively.  For  geometry  A  three  different  pitches  of  6.25mm,  12.7mm  and  25.4mm  were 
used.  The  diameter  of  the  holes  was  increased  in  steps  to  achieve  different  porosities. 
Porosity  is  defined  as: 


%  porosity  -  * [ ^/2 j  *  100 

2 

P 

where  d  is  the  diameter  and  p  is  the  pitch  of  the  holes. 

Due  to  the  constraints  in  machinability  of  Homalite  100  and  to  avoid  obscurity  of 
the  fringe  pattern  in  the  photographs,  the  porosity  was  limited  to  a  maximum  of  70%. 
Tables  1  and  2  list  the  pitch,  the  diameter  of  the  holes  and  the  corresponding  porosity 
which  were  used  with  geometry  A  experiments. 

For  geometry  B,  two  different  pitches  of  12.7mm  and  25.4mm  were  used.  Again, 
different  diameter  holes  were  used  for  each  of  these  pitches  to  achieve  different  porosities. 
Table  3  shows  the  pitch,  the  diameters  of  the  holes,  and  the  corresponding  porosity  which 
were  used  with  this  geometry. 

In  the  saturated  porous  media  experiments,  water  was  used  as  the  saturation  fluid. 
It  was  contained  within  the  pores  by  sandwiching  both  sides  of  the  photoelastic  models  by 
means  of  thin  plexiglass  sheets  [14].  The  sides  of  these  sandwiched  models  were  made 
leak  proof  by  pressing  them  with  pressure  binders.  Care  was  taken  to  minimize  the 
number  of  air  bubbles  so  as  to  achieve  a  high  degree  of  saturation  [15,16]. 

The  wave  propagation  phenomenon  due  to  explosive  loading  was  studied  using  the 
technique  of  dynamic  photoelasticity  and  high  speed  photography.  The  models  were  placed 
in  the  optical  bench  of  the  high  speed  multiple  spark  gap  camera.  The  camera  was 
triggered  at  a  prescribed  delay  time  after  igniting  the  explosive.  The  high  speed 
photographic  system  operates  as  a  series  of  extremely  short  duration  pulses  of  high 
intensity  light  and  provides  20  photoelastic  images  at  discrete  time  intervals  during  the 
dynamic  event.  These  photographs  of  the  wave  propagation  process  at  different  stages  of 

development  provided  the  necessary  data  to  obtain  the  velocity  and  attenuation  of  the 
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leading  from  of  the  longitudinal  wave. 


ANALYSIS 

The  sequence  of  20  photographs  obtained  from  the  high  speed  camera  were  analyzed 
to  obtain  the  compressional  wave  velocity  and  the  stress  wave  attenuadon  in  both  the 
unsaturated  and  the  saturated  porous  media.  The  compressional  wave  velocity  was 
obtained  by  plotring  the  instantaneous  position  of  the  wavefront  with  respect  to  time.  The 
gradient  of  this  line  gave  the  average  compressional  wave  velocity.  Figure  3  shows  a 
typical  position  versus  time  plot  The  portion  of  the  curve  between  X  and  Y  represents  the 
region  during  which  the  stress  wave  was  in  the  porous  zone.  The  portion  of  the  curve 
before  point  X  and  after  point  Y  represent  the  regions  before  the  stress  wave  entered  the 
porous  zone  and  after  it  emerged  from  it  respectively. 

The  fringe  patterns  were  further  analyzed  to  determine  the  propagation  and 
attenuation  characteristics  of  the  porous  media.  To  obtain  the  attenuation  characteristics,  a 
plot  of  the  maximum  normalized  fringe  order  around  the  hole  versus  the  hole  number 
(along  the  line  0=0)  was  obtained.  The  normalization  was  done  with  respect  to  the 
maximum  fringe  order  on  the  first  hole  directly  below  the  explosive  charge.  To  separate 
the  effects  of  material  and  geometric  attenuation,  an  expression  of  the  form 


N 

max 


for  the  maximum  fringe  order  was  sought,  where  y  is  the  distance  from  the  loading  source 
and  k,  the  coefficient  of  material  attenuation.  The  term  y'm  accounts  for  the  geometric 
attenuation  in  two  dimensional  space  (plate).  The  attenuation  coefficient  and  the  constant 
N0  were  obtained  by  plotting  NMly'l/2  vs.y  on  a  semilog  scale.  In  order  to  compare  the 
saturated  and  unsaturated  porous  media  with  the  unflawed  material  the  attenuation  curves 
for  all  three  were  plotted  together. 


The  first  series  of  experiments  was  conducted  with  geometry  A  and  involved  five 
models  with  pitch  equal  to  12.7mm.  The  diameter  of  the  holes  was  increased  in  steps  to 
achieve  porosity  ranging  from  zero  to  seventy  percent.  With  the  detonation  of  the 
explosive,  longitudinal  and  shear  waves  were  produced  which  travelled  with  velocities  of 
2140  m/s  and  1220  m/s  respectively  in  the  plate.  As  the  longitudinal  wave  passes  over  the 
first  row  of  the  array  of  holes  the  stresses  around  them  start  to  build  up.  Since  the 
diameter  of  the  holes  is  much  less  than  the  radius  of  curvature  of  the  approaching  wave 
front  the  fringes  tend  to  align  themselves  in  a  direction  which  is  normal  to  the  approaching 
wavefront  as  shown  in  Fig  4.  As  the  wavefront  impinges  upon  the  hole  boundary 
scattering  and  diffraction  takes  place  and  reflected  longitudinal  waves  ard  shear  waves  are 
produced.  The  reflected  waves  interact  with  the  propagating  wave  altering  the 
semicircular  na  ure  of  the  wavefront.  With  the  larger  diameter  holes,  the  ratio  X/d  (where 
X  is  the  wavelength  of  the  pulse  and  d  the  diameter  of  the  holes)  decreases  resulting  in  an 
increased  level  of  scattering  and  diffraction. 

Figure  5  shows  a  typical  set  of  photographs  obtained  during  the  wave  propagation 
experiments.  These  dynamic  isochromatic  fringes  show  the  development  of  the  wavefront 
as  the  wave  propagates  through  the  porous  medium.  With  continued  reflection  from  each 
row  of  holes  and  the  superposition  of  the  reflected  waves,  the  curvature  of  the  leading 
dilatational  wavefront  keeps  on  decreasing  as  the  wave  moves  further  into  the  porous  zone. 
This  phenomenon  increases  for  higher  porosities.  Figure  6  shows  the  location  and  the 
shape  of  the  wavefront  for  a  given  porosity. 

The  dynamic  photoelastic  photographs  were  analyzed  to  get  the  leading  wave 
velocity.  Figure  7  shows  the  variation  in  wave  velocity  as  a  function  of  porosity.  The 
dilatational  wave  velocity,  measured  in  a  direction  directly  below  the  charge  (along  the  line 
theta=0),  is  less  than  the  dilatational  wave  velocity  in  an  unflawed  half-plane.  The  wave 
velocity  decreases  by  almost  18%  as  the  porosity  is  increased  from  0  to  20%.  As  the 
porosity  is  increased  further,  the  drop  in  velocity  is  small  until  the  50%  porosity  mark  is 
passed,  after  which  wave  velocity  decreases  sharply  and  drops  down  to  almost  60%  of  its 
value  for  the  unflawed  half-plane. 

The  dilatational  wave  velocity  also  shows  directional  dependence  and  tends  to 
decrease  as  the  angle  theta  is  increased.  Figure  8  shows  that  the  velocity  at  any  angle 
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theta  depends  upon  the  diameter  of  the  holes  in  the  porous  zone.  For  small  diameter 
holes,  little  change  in  velocity  can  be  seen  as  the  angle  theta  is  increased  from  0-30%. 
But  as  the  diameter  of  the  holes  is  increased,  the  wave  velocity  drops  noticeably  and  this 
drop  increases  as  the  angle  theta  is  increased.  Around  0  equal  to  30°,  there  is  a  drop  in 
velocity  of  6%  for  the  9.6mm  diameter  holes,  10%  for  the  10.9  mm  diameter  hole  and 
16%  for  the  11.7mm  holes.  This  dependence  of  the  wave  velocity  in  the  porous  zone  on 
the  angle  theta  is  due  to  the  fact  that  the  wave  sees  a  different  kind  of  microstructure  for 
different  values  of  theta.  This  decrease  in  velocity  also  accounts  for  the  decrease  in  the 
curvature  of  the  wavefront  mentioned  earlier. 

Effect  of  Pitch 

A  series  of  experiments  was  conducted  to  study  the  effect  of  pitch,  or  the  spacing 

between  the  holes  on  the  wave  propagation  phenomenon.  Figure  9  shows  typical 

photographs  of  the  stress  wave  propagation  in  the  porous  medium  for  three  different  values 

of  pitch.  As  the  pitch  is  increased,  the  ratio  of  the  number  of  holes  to  that  of  the 

wavelength  of  the  stress  wave  decreases,  resulting  in  reduced  diffraction  and  scattering  of 
the  stress  wave.  The  influence  on  the  stress  around  the  holes  due  to  neighboring  holes  is 
also  reduced. 

Figure  10  shows  the  normalized  wave  velocity  vertically  below  the  explosive  charge 
as  a  function  of  porosity  for  holes  at  three  different  pitches.  For  all  three  pitches,  the 
wave  velocity  decreases  as  the  pov>sit'-  of  the  medium  is  increased.  However,  for  a  given 
porosity,  the  wave  velocity  increase.,  s  the  pitch  is  increased.  This  behavior  in  stress 
wave  velocity  is  also  observed  in  other  directions  (0  *  0). 

Effect  of  Geometry 

To  study  the  effect  of  microstructure  or  the  arrangement  of  the  pores  on  wave 
propagation  characteristics,  the  arrangement  of  the  holes  in  the  porous  media  was  changed 
from  geometry  A  (cubic  arrangement)  to  geometry  B  (hexagonal  arrangement).  Figure  11 
shows  typical  isochromatic  fringes  obtained  during  wave  propagation  in  geometry  B. 
These  fringes  show  more  scattering  and  diffraction  than  obtained  with  geometry  A.  Figure 
12  shows  the  wave  velocity  as  a  function  of  porosity  vertically  below  the  explosive  charge 
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for  geometry  B.  The  velocity  in  geometry  B  decreases  with  porosity  like  it  does  in 
geometry  A.  The  velocity  increases  as  the  pitch  goes  up  for  a  given  porosity.  Although 
the  velocity  vertically  below  this  explosive  charge  shows  a  similar  trend  in  both 
geometries,  the  absolute  value  of  velocity  in  geometry  B  is  much  lower  than  in  geometry 
A.  This  comparison  is  shown  in  Fig.  13. 

For  geometry  B,  the  stress  wave  velocity  also  shows  directional  dependence  but 
here,  unlike  geometry  A,  the  velocity  increases  as  the  angle  is  increased  from  0-30°. 

Figure  14  shows  the  normalized  wave  velocity  as  a  function  of  angle  ©  for  a  given  pitch 
and  varying  porosities  in  geometry  B.  The  velocity  increases  with  0  and,  for  0  equal  to 
30°,  'his  velocity  almost  approaches  the  velocity  for  geometry  A  at  0=0°.  A  similar  trend, 
but  in  the  opposite  direction,  can  be  seen  for  geometry  A  where  the  stress  wave  velocity 
decreased  with  angle  0  and  for  0=30°  approached  the  wave  velocity  in  geometry  B  at 
0=0°.  This  behavior  can  be  explained  by  looking  closely  at  the  microstructure  the  wave 
sees  as  the  angle  0  is  increased  in  both  geometries.  Figure  15  shows  a  schematic  of  the 
microstructure  the  wave  sees  at  an  angle  of  45°  in  geometries  A  and  B.  It  is  interesting  to 
note  that,  at  45°,  the  pore  arrangement  in  geometry  B  appears  like  geometry  A  at  0°  and 
vice-versa.  Changing  the  angle  further  from  45*  to  90“  causes  the  microstructure  to  go 
back  to  the  one  observed  for  0=0°. 

Effect  of  Moisture  in  the  Pores 

A  series  of  experiments  was  conducted  with  geometry  A  to  study  the  effect  of 
moisture  in  the  pores  on  the  wave  propagation  phenomenon.  Figure  16  shows  typical 
isochromatic  fringes  obtained  when  the  stress  wave  propagates  in  porous  media.  These 
photographs  were  analyzed  to  study  both  the  wave  velocity  as  well  as  the  stress  wave 
attenuation.  Figure  17  shows  the  normalized  wave  velocity  as  a  function  of  the  porosity 
for  three  different  pitch  arrangements.  The  trend  in  velocity  is  very  similar  to  the  one 
obtained  for  unsaturated  media.  However  the  velocities  in  the  saturated  media  are  much 
lower  than  those  for  unsaturated  media  as  shown  in  Fig.  18.  For  small  porosities  (r\  < 
20%)  the  amount  of  moisture  per  unit  volume  is  very  small  and  therefore  there  is  no 
difference  in  velocities.  As  the  porosity  is  increased  beyond  20%  the  velocity  shows  a 
much  more  rapid  drop  in  a  saturated  media  than  an  unsaturated  one.  This  can  be 
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explained  to  some  extent  by  the  fact  that  the  composite  density  for  the  porous  media  with 
water  is  higher  than  that  of  a  porous  media  with  air  in  the  holes.  These  results  are  in 
agreement  with  those  obtained  by  others  [1, 2,7,9]  for  saturated  soils. 

Attenuation  of  the  Stress  Waves 

The  compressive  wave  front  attenuates  rapidly  as  the  stress  wave  propagates  through 
the  porous  media.  As  in  the  case  of  the  stress  wave  velocity,  the  attenuation  of  the  stress 
wave  depends  strongly  on  the  microstructure  of  the  porous  media.  Larger  the  size  of  the 
holes,  the  higher  is  the  porosity  and  greater  is  the  attenuation  of  the  stress  wave. 
Moreover,  the  wave  shows  a  higher  attenuation  as  the  pitch  or  the  spacing  between  the 
holes  is  reduced.  The  geometry  or  the  configuration  of  the  holes  also  effects  m  the 
attenuation  characteristics  and  there  is  higher  attenuation  for  Geometry  "B"  than  for 
Geometry  "A".  This  is  to  be  expected  since  the  stress  wave  undergoes  much  more 
scattering  and  diffraction  when  it  propagates  through  the  porous  media  with  Geometry  "B" 
than  it  does  in  Geometry  "A". 

Figure  19  shows  the  isochramatic  fringe  pattern  showing  the  circumferential  fringe 
distribution  around  the  periphery  of  the  holes  for  two  different  hole  diameters.  To  obtain 
the  attenuation  characteristics,  a  plot  of  the  maximum  normalized  fringe  order  vs.  the  hole 
number  (along  the  line  theta=0)  was  obtained.  Figure  20  shows  the  plot  of  the  attenuation 
curves  for  both  the  unsaturated  and  saturated  porous  media  for  Geometry  "A”  and  pitch  = 
12.7mm.  The  attenuation  curves  are  drawn  for  porosity  of  11%  and  20%.  Figure  21 
shows  the  attenuation  curves  for  the  case  of  unsaturated  porous  media  with  Geometry  B 
and  porosity  equal  to  20%  and  46%,  respectively.  For  the  sake  of  comparison,  the 
geometric  attenuation  in  an  unflawed  halfplane  is  also  drawn  on  the  same  plots. 

To  characterize  the  attenuation  of  the  stress  wave  in  the  porous  media  and  to 
separate  the  geometric  damping  from  the  material  damping,  was  plotted  as  a 

function  of  y  on  a  semilog  scale  as  shown  in  Figs  23  &  24.  These  plots  were  analyzed  to 
calculate  the  attenuation  coefficients  which  are  shown  in  Table  4.  These  values  give  us 
an  idea  of  the  attenuation  characteristics  of  the  porous  material.  Higher  values  of  K 
indicate  higher  attenuation  of  the  stress  wave  in  that  material. 
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Dynamic  photoelasticity  along  with  high  speed  photography  were  used  to  study  wave 
propagation  in  porous  media.  Plate  specimens  were  fabricated  from  Homalite-100  and  an 
array  of  holes,  having  a  specific  geometry,  were  machined  in  them  to  simulate  the  porous 
media.  The  specimens  were  loaded  explosively  on  one  edge  and  the  resulting  isochromatic 
fringe  patterns  were  analyzed  to  obtain  the  compressional  wave  velocity  and  the  stress 
wave  attenuation  in  the  porous  material. 

The  tress  wave  speed  dropped  with  increasing  porosity.  For  a  given  porosity,  the 
wave  speed  is  lower  for  smaller  pitch  in  both  the  geometries.  The  wave-velocity  also 
showed  strong  dependence  on  the  micro  structure  of  the  porous  media.  Geometry  B  showed 
a  larger  drop  in  velocity  than  Geometry  A  for  the  same  porosities.  Due  to  the  point 
loading  and  the  geometry  of  the  porous  models,  the  stress  wave  exhibits  directional 
dependence  as  it  sees  a  different  kind  of  microstructure  as  a  function  of  the  angle  theta. 

The  stress  wave  shows  appreciable  attenuation  in  the  porous  media.  The  attenuation 
observed  was  separated  into  its  geometric  and  material  components.  The  coefficients  for 
the  material  attenuation  were  calculated  for  different  porous  medias  and  are  summarized  in 
Table  4.  From  the  table  it  can  be  seen  that  the  attenuation  constants  are  higher  for  the 
case  of  porous  media  with  a  higher  porosity  and  they  decrease  for  the  porous  media  with  a 
larger  pitch  but  having  the  same  porosity.  Moreover  the  attenuation  constants  are  higher 
for  the  configuration  of  the  holes  in  geometry  B  than  in  geometry  A. 

The  wave  propagation  characteristics  also  depend  on  whether  the  porous  medium 
is  nearly  saturated  or  unsaturated.  Nearly  saturated  porous  media  show  a  larger  drop  in 
wave  velocity  than  the  unsaturated  media.  Moreover  the  stress  wave  attenuation  is  much 
higher  in  a  nearly  saturated  medium  than  an  unsaturated  one. 
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travels  through  saturated  and  unsaturated  porous 
media  of  geometry  A. 
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DYNAMIC  PHOTOELASTIC 
INVESTIGATION  OF  WAVE 
PROPAGATION  AND  ENERGY 
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H.  P.  ROSSMANITH  Institute  of  Mechanics,  Technical  University,  Vienna,  Austria 

This  paper  deals  with  the  dynamic  contact  of  an  explosively  excited  disc  with  another  disc  or  a  half-plane.  Dynamic 
photoelastic  recordings  show  the  development  of  the  time-dependent  contact  area  and  the  formation  of  the  highly 
complex  diffraction  pattern. 


1  INTRODUCTION 

The  dynamic  contact  between  elastic  bodies  and  the 
transfer  of  momentum  and  energy  across  the  contact 
area  is  fundamental  to  studies  of  dynamic  load  transfer 
in  granular  or  particulate  media.  When  a  body  made  up 
of  a  large  number  of  cohesionless  particles  is  impacted 
upon  with  a  punch  or  an  explosive,  transient  dynamic 
load  transfer  paths  develop  which  differ  from  those 
established  during  static  loading.  This  difference  occurs 
because,  in  dynamic  loading,  inertia  effects  play  an 
important  role.  Most  of  the  research  to  date  in  contact 
mechanics  pertains  either  to  static  studies  or  when 
dynamic  impact  is  considered;  first  order  discretized 
lump  mass  models  are  employed.  The  dynamics  and 
mechanics  of  physical  impact  between  solid  bodies  up  to 
the  late  1960s  is  covered  in  references  (lH4)t-  Generally 
in  these  studies  the  impacting  body  is  considered  rigid 
and  stress  wave  effects  in  the  indentor  are  not  taken  into 
account.  Also,  in  collision  problems,  the  bodies  are 
treated  as  particles  and,  again,  wave  propagation  pro¬ 
cesses  occurring  within  the  body  are  not  taken  into 
account.  Recently  studies  by  Comninou  et  al.  (5X6)  have 
focussed  on  the  transmission  of  wave  motion  between 
two  solids.  In  these  theoretical  studies  plane  waves  are 
considered  passing  from  one  half  plane  to  another. 
Several  initial  conditions  are  looked  at,  including  an 
initial  gap  and  also  friction  at  the  interface.  A  recent 
experimental  study  of  dynamic  load  transfer  in  a  granu¬ 
lar  media  was  done  by  Rossmanith  and  Shukla  (7)  where 
load  transfer  paths  in  both  systematic  and  random 
aggregates  of  disc  assemblies  were  recorded  using 
dynamic  photoelasticity. 

This  paper  addresses  the  problem  of  formation  of 
dynamic  contacts  between  two  solids.  The  experimental 
technique  of  dynamic  photoelasticity  is  utilized  to  visual¬ 
ize  the  formation  of  contact  and  wave  scattering  between 
the  bodies. 

The  MS.  of  this  paper  was  received  at  the  Institution  on  9  September  !98S  and 
accepted  for  publication  on  13  May  1986 
t  References  are  given  in  the  Appendix 
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2  EXPERIMENTAL  PROCEDURE 
The  experimental  method  utilized  in  this  study  was 
dynamic  photoelasticity  along  with  high  speed  pho¬ 
tography.  Dynamic  photoelasticity  provides  global  field 
data  during  the  dynamic  contact,  such  as  the  change  of 
the  state  of  stress  in  the  vicinity  of  the  contact  region  and 
in  the  bulk  of  the  contacting  solids,  as  a  function  of  time 
and  space. 

In  order  to  understand  the  physics  of  dynamic  contact 
between  two  elastic  solids  a  series  of  four  controlled 
experiments  was  performed,  where  an  explosively  loaded 
disc  was  in  contact  with  another  disc  or  a  half-plane  and 
the  explosive  excitation  was  located  either  at  the  centre 
of  the  disc  or  at  the  rim.  Figure  1  shows  the  photoelastic 
model  and  explosive  arrangement  and  Table  1  lists  the 
characteristics  of  the  dynamic  loading. 

The  circular  discs  of  radius  R  =  76.2  mm  and  the  half¬ 
plane  were  fabricated  from  a  6.35  mm  thick  sheet  of 
Homalite  100  which  is  a  commercially  available  clear 
transparent  polyester.  This  photoelastic  material 
becomes  temporarily  birefringent  when  subjected  to  a 
state  of  stress  and  gives  rise  to  optical  interference  fringes 
when  viewed  in  a  circular  polariscope.  These  fringes  are 
known  as  isochromatics  and  represent  lines  along  which 
the  difference  of  the  principal  normal  stresses  is  constant 

Dynamic  loading  was  achieved  by  means  of  small 
charges  of  150  mg  PETN  (experiments  Nos  1  and  2)  or  of 
100  mg  lead  azide  (PbN6)  (experiments  Nos  3  and  4).  The 
explosive  excitation  source  is  located  at  the  centre.  A,  of 
the  disc  in  experiments  Nos  1  and  2  and  at  the  top,  B,  in 
experiments  Nos  3  and  4.  In  the  experiments  with  a 


Table  1.  List  of  experiments  of  dynamic  contact 


Experiment 

number 

Configuration 

Location  of 
charge 

Charge 

No.  1 

Disc/half-plane 

Centre 

1 50  mg  PETN 

No.  2 

Disc/disc 

Centre 

150  mg  PETN 

No.  3 

Disc/disc 

Top 

100  mg  PbN6 

No.  4 

Disc/half-plane 

Top 

100  mg  PbN6 

HI 


centre  charge,  a  specially  developed  pressure  contain¬ 
ment  device  was  utilized  to  prevent  the  field  of  view  from 
being  obscured  by  the  rapidly  expanding  cloud  of  com¬ 
bustion  gas  products. 

For  photoelastic  fringe  pattern  recording,  a  multiple 
spark  gap  camera  of  the  Cranz-Schardin  type  was  used. 
The  system  and  its  use  in  dynamic  photoelasticity  have 
been  documented  frequently  in  the  literature,  see,  e.g.,  (9). 
The  experimental  outcome  is  a  record  of  a  time- 
controlled  sequence  of  20  images  of  the  spatial  configu¬ 
ration  of  the  dynamic  process  at  discrete  times.  Thus,  a 
sequence  of  20  ultra-short  flash  photographs  are 
obtained  that  show  individual  stages  of  the  evolution  of 
the  wave  interaction  process.  In  these  wave  propagation 
experiments  the  camera  was  operated  at  200000  frames 
per  second  and  provided  an  observation  period  of  100  ns 
with  an  additional  initial  deiay  of  25  ^s. 

3  THEORETICAL  CONSIDERATIONS 

3. 1  Size  of  dynamic  contact  area 

Consider  two  discs  cf  radii  R,  (i  =  1,  2)  touching  at  a 
point  PC(RV  0)  as  shown  in  Fig.  1.  Let  the  centres  of  the 
discs  be  fixed  at  Af,( 0,0)  and  Mj(R,  +  R2,0).  Upon 
detonation  of  a  cylindrical  explosive  charge  at  the  centre 
of  disc  No.  1  a  circular  crested  detonation  pulse  radiates 
from  A#,.  The  displacements,  ur,  u,,  and  the  stresses,  cr„, 
Oft  vary  according  to  Bessel  functions.  For  ‘large’  values 
of  r,  the  asymptotic  expansions  yield  a  plane  wave 
approach.  ‘Large  values’  of  r  are  attained  within  4  to  5 
zeroes  of  the  Bessel  function,  and  in  the  experiments  this 
distance  r  corresponds  to  a  few  plate  thicknesses  b  (r  % 
86).  Furthermore,  any  explosively  induced  pulse  dis¬ 
turbance  can  be  considered  the  result  of  the  super¬ 
position  of  a  set  of  harmonic  waves. 

The  data  reduction  procedures  for  circular  cylindrical 


pressure  waves  as  developed  in  (II)  render  relations 
between  the  strains,  er,  e„,  and  stresses,  <rr,  <x#,  and  the 
isochromatic  fringe  order  distribution.  N.  In  2D  dynamic 
photoelastic  wave  propagation  experiments  an  r~' 
attenuation  as  associated  with  an  expanding  spherical 
wave  is  commonly  observed  which  is  in  marked  contrast 
to  the  r~ 1/2  attenuation  of  the  amplitude  for  cylindrical 
waves  as  predicted  by  theory  (12X13).  For  wave  positions 
far  from  the  detonation  centre  a  linear  relationship 
between  radial  stress,  <xr,  and  radial  particle  velocity,  u, 
(14),  at  =  pc,t\  may  be  assumed  to  hold  throughout  the 
pulse.  This  ‘far-field’  approximation  corresponds  to  a 
local  replacement  of  the  detonation  wave  by  a  plane 
wave  with  r~l  attenuation  and  is  supported  by  experi¬ 
mental  findings. 

For  contact  investigations,  the  radial  expansion  of  the 
edge  of  the  explosively  excited  disc  is  of  interest.  Free 
radial  expansion  can  occur  along  the  circumference 
except  in  the  region  of  contact.  In  Fig.  1,  as  disc  No.  1 
expands,  the  contact  area  increases  in  time.  The  dynamic 
contact  problem  is  further  complicated  due  to  pulse  dif¬ 
fraction  at  the  edge  of  contact  and  transmission  across 
the  contact  due  to  time  dependent  boundary  conditions. 
The  situation  is  illustrated  in  Fig.  2,  where  two  different 
stages  of  contact  area  development  during  the  passage  of 
a  half  sine  pulse  are  shown.  In  Fig.  2  the  distance  A' A"  of 
the  intersection  points  of  the  two  circles  with  radii  R, 
4-  ur  and  R2  is  an  upper  limit  for  the  extension  of 
contact.  For  two  discs  the  distance  A'A"  is  given  by 

A' A"  =  2a(f)  =  IRyJ^l  +  70  -  £jj  (1) 
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Fig.  2.  Dynamic  contact  of  two  discs  showing  contact  area  develop¬ 
ment 


is  shown  in  Fig.  3.  Time-dependent  wave  diffraction 
about  the  moving  contact  edge  during  increasing  and 
receding  contact  (contact  cycles)  poses  a  challenging 
theoretical  problem.  An  exact  quantitative  analytical 
treatment  of  this  problem  is  not  feasible  to  date.  Numeri¬ 
cal  work  for  simulation  of  dynamic  contact  behavior  is  in 
progress. 

4  EXPERIMENTAL  RESULTS  AND  DISCUSSION 
The  results  of  the  dynamic  photoelastic  experiments  will 
be  discussed  here.  For  convenience  the  experiments  with 
centre  charge  and  the  top  charge  will  be  discussed 
separately. 

4.1  Dynamic  contact  experiments  with  centre  charge 

The  first  experiment  refers  to  the  disc-half  plane 
geometry.  An  explosive  charge  of  150  mg  of  PETN  was 
detonated  at  the  centre  of  the  disc.  A  sequence  of  three 
photographs  obtained  during  the  experiment  are  shown 
in  Fig.  4.  Upon  detonation  of  the  explosive  charge  a 


which  reduces  to  e  =  1  for  the  disc  -  half  plane  experi¬ 
ments.  Results  for  Homalite  100  with  c,  -  2120  m,  sec,  a 
half  sine  pulse  length  of  A  *  =  12.5  mm  and,  of  peak 
fringe  order  ,Vm„  =  30  for  the  disc-disc  arrangement  (R  ( 
=  R2  =  R  =  76.2  mm)  and  the  half  plane-disc  arrange¬ 
ment  (R |  =  76.2  mm.  R2  =  x).  yield  contact  times  of 
1 2  //s  and  maximum  upper  bound  for  the  contact  area  of 
a  R  $  0.007  (disc-half  planel  and  iiR$  0.005  (disc- 
disc) 

Time  change  of  size  of  contact  area  depends  primarily 
on  the  slope  of  the  incident  wave.  Once  contact  is  estab¬ 
lished  for  an  elemental  section  of  the  circumference, 
energy  is  transferred  to  the  receptor  disc  and  the  ratios  of 
particle  speeds  of  the  respective  waves  involved  in  reflec¬ 
tion  and  transmission  across  the  contact  area  are  the 
same  as  the  ratios  between  the  wave  amplitudes.  A 
typical  wave  front  construction  for  the  dynamic  contact 


Fig.  3.  Diffraction  and  reflection  of  3n  incident  elastic  P  wave  during 

formation  of  elastic  contact  - reflected  waves.  - - transmitted 

waves 


Frame  8.  t  =  56 


Frame  10.  t  =  65  ys 


Frame  1 8.  i  =  1 1 3 


Fig  4  Centre  charge  induced  dynamic  contact  showing  P  wave  inter¬ 
action 
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cylindrical  longitudinal  (compressive)  wave,  P,.  emerges 
and  radiates  from  the  centre  of  the  disc.  Due  to  cracking 
and  breakdown  of  the  borehole  the  rotational  symmetry 
is  disturbed  and  this  gives  rise  to  a  P  wave  along  with  a 
random  shear  wave  which  quickly  outdistances  the 
emerging  cracks.  Hence.  P  and  S  waves  are  emitted  from 
the  detonation  site  and  interact  with  the  contact  zone. 

The  first  photograph,  frame  No.  8.  taken  at  t  =  56  /is 
after  detonation,  shows  the  dynamic  situation  at  a  time 
shortly  after  initial  contact  is  established.  The  expanding 
compressive  incident  P,  wave  has  been  reflected  along 
the  free  sections  of  the  boundary  to  yield  a  tensile  P,P, 
wave.  Part  of  the  energy  and  momentum  of  the  wave  has 
been  transferred  across  the  contact.  The  exact  location  of 
the  front  of  diffracted  transmitted  P:  P,  wave  is  difficult 
to  locate  (121;  however,  it  is  a  common  practice  to  iden¬ 
tify  the  wave  front  with  a  line  drawn  just  ahead  of  the 
V  =  1  2  fringe.  Although  the  contact  area  is  obscured  by 
a  pseudocaustic,  an  estimate  of  the  contact  length  was 
made.  At  time  r  =  56  /is  the  contact  length  was  approx¬ 
imately  1 5  mm. 

In  frame  No.  10  of  Fig.  4.  taken  at  t  =  65  /is,  the 
reflected  and  transmitted  longitudinal  and  shear  pulses 
have  separated  and  the  resulting  fringe  pattern  in  the  half 
plane  is  similar  to  the  pattern  that  would  be  generated 
upon  detonation  of  a  concentrated  charge  at  the  centre 
of  the  contact  zone.  With  increasing  time  after  contact 
the  fringe  pattern  in  the  half  plane  unfolds  into  a  longitu¬ 
dinal  wave  P,P,.  a  shear  wave  S;P,,  a  von-Schmidt  or 
head  wave.  V,  P,.  and  a  Rayleigh  surface-wave,  R,  P,. 
The  corresponding  waves  in  the  disc  are  labeled  P,P,, 
S,P,.  V,P,,and  R , P , .  Rayleigh-wave-induced  cracking 
in  the  disc  at  a  position  adjacent  to  the  end  of  the  contact 
zone  was  observed  in  this  experiment. 

Progressive  dynamic  contact  formation  due  to  shear 
wave  interaction  with  the  contact  is  shown  in  the  third 
photograph  of  Fig.  4.  frame  No.  18  taken  at  f  =  113 
/isec.  The  front  of  the  S;  S,  wave  is  located  43  mm  below 
the  half  plane  boundary  and  this  matches  well  with  the 
theoretical  predictions.  The  R2  P,  Rayleigh  waves  have 
completely  separated  from  the  other  bulk  waves  and  pro¬ 
pagate  with  speed  c„  =  110  m  sec  along  the  free  bound¬ 
ary.  The  increased  number  of  fringes  at  the  contact  zone 
is  indicative  of  increased  contact  stresses. 


Frame  7,  r  =  P,4 


Fig  5  Top  charge  induced  dynamic  contacl  showing  P  wave  inter¬ 
action 


Frame  12.  r=  127  jas 


Frame  15.  r=  155  ps 


Fig.  6.  Top  charge  induced  dynamic  conlacl  showing  S  wave  inter 
action 

In  a  second  experiment  a  disc  assembly  was  used  with 
the  same  centre  charge  The  results  were  qualitatively 
similar  to  the  ones  obtained  in  experiment  No.  1 . 

4.2  Dynamic  contact  experiments  with  top  charge 

The  third  experiment  was  done  with  a  disc-disc 
assembly.  A  charge  of  100  mg  of  lead  azide  was  exploded 
on  the  top  edge  of  the  first  disc.  This  explosion  generates 
P  and  S  waves  which  propagate  in  the  bodv  of  the  disc 
and  Rayleigh  surface  waves  that  travel  along  the  circum¬ 
ference  of  the  disc. 

The  P  wave  which  travels  with  a  velocitv  oft',  =  21 20 
m/sec  is  the  first  wave  to  interact  with  the  contact  This 
interaction  is  shown  in  Fig  5. 

Frame  No.  7,  taken  at  r  —  84  /is  after  detonation, 
shows  the  establishment  of  contact  where  part  of  ihe 
wave  energy  is  already  transferred  across  the  contact. 
The  contact  area  is  proportional  to  the  size  of  the  pseudo 
caustic  that  forms  at  the  contact  zone  (15).  The  contact 
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Frame  13.  /  =  20 1  ps 


Fig.  7.  Top  charge  induced  dynamic  contact  showing  R  wave  inter¬ 
action 


length  in  this  frame  is  estimated  to  be  approximately 
1.5  mm.  Diffraction  and  transmission  of  the  incident  P 
wave  generates  P  and  S  waves  in  both  discs  and  also 
gives  rise  to  R  waves  that  travel  along  the  outer  bound¬ 
aries  of  the  discs. 

The  sequence  of  photographic  recordings.  Fig.  6, 
shows  the  interaction  of  the  incident  S  wave  with  the 
contact.  Frame  No.  12,  taken  at  t  =  127  /is,  depicts  the 
early  stage  of  S  wave  diffraction  about,  and  S  wave  trans¬ 
mission  across,  the  contact  surface.  Due  to  the  shearing 
nature  of  the  wave  the  fringe  pattern  in  the  vicinity  of  the 
contact  region  is  unsymmetrical.  Moreover,  the  experi¬ 
mental  recordings  reveal  that  this  fringe  pattern  oscil¬ 
lates  about  the  normal  to  the  contact  plane.  This  effect 
can  be  seen  in  later  frames.  In  frame  No.  15  at  t  =  155  /gs 
the  smallest  visible  fringe  below  the  contact  in  the  recep¬ 
tor  disc  is  of  order  /V  =  6  as  opposed  to  a  maximum 
visible  fringe  order  of  N  =  3  in  frame  No.  12.  This 
increase  is  indicative  of  severe  contact  stresses.  This  is 
also  evident  from  the  increase  of  the  size  of  the  pseudo¬ 
caustics  at  the  edges  of  the  contact  area.  In  frame  No.  15 
the  contact  length  is  about  2  mm.  Frame  No.  16  at 
t  =  175  #js  shows  the  propagation  of  the  S2  S,  wave  and 
the  R2S,  Rayleigh  wave  generated  by  shear  wave  inter¬ 
action.  The  Rayleigh  wave  interaction  with  the  contact 
could  not  be  studied  in  this  experiment  as  the  time  dura¬ 
tion  of  the  experiment  was  not  long  enough. 
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The  second  experiment  was  performed  with  the  disc  - 
half  plane  assembly.  The  timings  were  adjusted  so  that  R 
wave  interaction  with  the  contact  could  be  observed.  The 
interaction  of  the  P  and  S  wave  with  the  contact  were 
similar  to  the  previous  experiment.  The  R  wave  inter¬ 
action  is  shown  in  the  sequence  of  photographs  in  Fig.  7. 
Frame  No.  13,  taken  at  t  =  210  /js,  shows  the  R  waves  as 
they  approach  the  contact  region.  Since  the  leading  part 
of  the  R  wave  exerts  circumferential  tension,  the  contact 
area  will  recede  due  to  radial  contraction.  This  is  shown 
in  frame  No.  14,  taken  at  t  =  215  jjs,  where  the  leading 
peaks  of  the  R  waves  overlap  and  the  contact  has  been 
lost.  In  frame  No.  15,  taken  at  t  =  232  /is,  ihe  compres¬ 
sive  trailing  parts  of  the  R  waves  superimpose  and 
contact  is  re-established.  Part  of  the  energy  of  the  R 
waves  is  transferred  into  the  half  plane,  as  can  be  seen 
from  the  fringe  development  in  frame  No.  15. 

5  CONCLUSIONS 

Dynamic  photoelastic  investigation  of  the  interaction  of 
elastic  stress  waves  with  contacts  in  solids  provides  full- 
field  information  of  diffraction,  reflection,  and  transmis¬ 
sion  processes  as  a  function  of  time.  Qualitative 
evaluation  of  dynamic  isochromatic  fringe  patterns 
reveals  that  the  individual  wave  types  interact  with  the 
contact  region  in  a  distinctive  way. 

Incident  P  waves  give  rise  to  relatively  strong  trans¬ 
mitted  and  reflected  P,  S,  and  R  waves.  The  intensity  of 
the  r  lected  R  waves  is  high  enough  to  initiate  cracking 
on  either  side  of  the  contact  zone  in  the  donor  disc.  This 
occurs  due  to  the  joint  action  of  large  tensile  stresses  at 
the  surface  in  the  leading  part  of  the  R  wave  and  high 
friction  in  the  contact  area. 

During  S  wave  interaction  with  the  contact  compara¬ 
tively  little  energy  and  momentum  is  transferred  across 
the  contact  area.  The  interface  shearing  stress  and  shear 
deformations  reverse  sign  during  S  wave  interaction. 

Very  little  energy  and  momentum  are  transferred  into 
the  receptor  solid  during  the  R  wave  interactions  because 
of  receding  contact  during  the  passage  and  superposition 
of  the  leading  R  pulse.  Although  contact  is  re-established 
during  the  passage  of  the  trailing  compressive  R  pulse, 
the  normal  component  of  displacement  is  very  small  and, 
consequently,  little  energy  is  transferred  across  the 
contact. 

Further  work  on  the  quantitative  evaluation  of  experi¬ 
mental  fringe  pattern  recordings  on  the  basis  of  the 
approximate  wave  analysis  developed  in  the  paper  is  in 
progress.  Attention  is  focused  on  quantities  of  interest, 
such  as  distribution  of  contact  stresses  as  a  function  of 
time  and  space,  and  amount  of  energy  and  momentum 
transferred  across  the  contact  area.  The  experimentally 
recorded  patterns  will  serve  as  input  for  numerical  simu¬ 
lation  of  dynamic  contact  behaviour.  A  hybrid 
experimental-numerical  analysis  of  contact  stress  fringe 
pattern  for  a  quasi  static  problem  where  the  wave  length 
is  larger  than  the  particle  size  has  already  been  developed 
(16). 
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Abstract 


An  experimental  investigation  is  conducted  to  study  dynamic  load  transfer  in  granular 
rock  media.  The  granular  media  is  modelled  as  a  one-dimensional  chain  of  disks  fabricated 
from  four  different  types  of  white  Vermont  marble.  The  study  mainly  focusses  on  the  effect 
of  the  microstructure  on  transient  pulse  propagation.  The  transient  pulse  is  generated  by 
exploding  a  small  charge  of  lead  azide  on  top  of  the  disk  assembly.  During  wave  propagation 
dynamic  contact  strains  are  recorded  using  strain  gages.  This  information  is  used  to 
calculate  wave  velocity  and  attenuation  as  a  function  of  the  cumulative  damage  in  the  disk 
assembly. 

Introduction 

The  response  of  particulate  materials  to  impulsive  loading,  has  been  of  substantial 
interest  to  both  engineers  and  researchers.  The  stress  wave  propagation  in  such  a  media, 
due  to  explosive  loading,  depends  primarily  upon  the  load  transfer  process  by  which  the 
mechanical  signals  are  transmitted.  This  phenomenon  is  related  to  the  properties,  geometry 
as  well  as  the  microstructural  arrangement  of  the  particles  in  the  media.  As  the  stress  wave 
passes  through  the  system  rearrangement  of  the  particles  takes  place.  Further,  depending 
upon  the  amplitude  of  the  stress  wave,  damage  can  occur  in  the  particles.  Thus  the 

transmitted  stress-wave  carries  along  with  it  the  information  regarding  the  microstructural 
re-arrangenents  of  the  grains.  The  wave  propagation  characteristics  of  the  resulting  red : j 

is  different  as  compared  to  that  of  the  virgin  material.  Moreover,  it  has  been  observed  that 
for  real  earth  materials  like  sand  and  rock,  this  load  transfer  phenomenon  is  also  a  function 
of  the  state  and  the  history  of  the  applied  load. 

The  concept  of  modelling  granular  media  as  an  array  of  elastic  particles  led  to  the 
initial  attempts  at  predicting  wave  propagation  through  such  media.  Early  work  by  Iida 
[1,2],  Hughes  and  Cross  [  3  ]  ,  Hughes  and  Kelly  [  4  ]  ,  Gassman  and  Brandt  [  5  ]  employed  a  norma  1 
granular  contact  force  concept.  An  excellent  review  by  Deresiewicz  [6]  summarizes  both  the 
static  and  dynamic  studies  prior  to  1958.  Experimental  studies  of  the  load  transfer  in 
granular  media  have  been  previously  conducted  by  Drescher  and  DeJosselin  [7]  who  simulated 
granular  media  by  using  assemblies  of  circular  disks.  They  studied  the  static  load  transfer 
through  the  assembly  of  the  disks  using  photomechanics. 

Sossmanith  and  Shukla  [8]  extended  this  idea  to  the  dynamic  case  with  the  help  of  high 
speed  photography.  Additional  dynamic  work  was  also  carried  out  by  Shukla  and  Damania  [9' 
and  Shukla  and  Zhu  [10]. 


Ii 


The  present  paper  focuses  attention  on  the  wave  propagation  and  load  transfer  in  a 
single  one  dimensional  chain  of  disk  assembly.  The  disks  were  fabricated  from  four  different 
grades  of  white  Vermont  marble  each  having  a  different  microstructure.  The  stress  wave  is 
generated  by  exploding  a  small  amount  of  Lead  Azide  on  top  of  the  disk  assembly.  The 
dependence  of  the  stress  wave  velocity  on  the  microstructure  of  various  grades  of  marble 
rocks  is  studied.  Attempt  has  also  been  made  to  investigate  and  quantify  the  damage  that 
occurs  in  these  one  dimensional  aggregate  of  disks  under  repeated  stress  wave  loading. 
Tensile  splitting  tests  are  performed  to  estimate  the  residual  strength  of  the  damaged  disks. 
Photomicrographs  are  taken  at  various  stages  of  the  shock  wave  loading  process.  Moreover, 
the  stress-wave  velocity  as  well  as  wave  attenuation  characteristics  are  obtained  as  a 
function  of  the  accumulated  damage  in  the  grains. 

The  results  indicate  that  a  correlation  exists  between  the  grain  size  and  the  stress- 
wave  velocity.  It  is  seen  that  the  average  stress-wave  velocity  is  higher  for  rock  disks 
with  larger  gram  size.  Rapid  attenuation  in  the  stress  wave  peak  load  occurs  in  the  first 
few  disks  and  then  attains  a  steady  rate  of  decay.  The  residual  tensile  strength  of  the 
disks  first  decreases  with  the  applied  load  and  then  approaches  a  limiting  value  before 
eventually  fracturing.  The  stress-wave  velocity  rises  considerably  and  reaches  a  terminal 
velocity  upon  repeated  shock  wave  loading.  The  photomicrographs  taken  at  various  stages  of 
the  shock  wave  loadings  indicate  that  large  scale  crumbling  and  pitting  takes  place  near  the 
contact  zone.  This  leads  to  microcracking  and  eventual  failure  of  the  disks. 


Experimental  Procedure  and  Analysis 


A  series  of  experiments  are  conducted  using  a  single  chain  assembly  of  disks  as  shown 
in  Fig.  1.  These  disks  were  fabricated  from  four  different  cores  of  Vermont  marble  having 
different  microstructures  as  shown  in  Fig.  2.  The  disks  were  1.25  inches  (31.75mm)  in 
diameter  and  0.5  inches  (12.7mm)  in  thickness.  The  assembly  of  the  disks  was  explosively 
loaded  using  15mg  of  Lead  Azide  in  a  specially  designed  steel  charge  holder.  The  resulting 
dynamic  load  transfer  phenomenon  in  the  disks  was  studied  using  electrical  resistance  strain 
gages.  In  the  experiments  the  wave  length  of  the  explosive  loading  pulse  was  sufficiently 
large  as  compared  to  the  disk  diameters,  thus  resulting  in  a  quasi-static  type  of  loading 
around  the  contact  zone.  Hertz  contact  theory  along  with  the  experimentally  obtained  strains 
were  used  to  obtain  the  contact  stresses,  strains  and  the  dynamic  loads. 


From  the  Hertz  contact  stress  theory,  the  stress  field  equations  around  the  contact 
region  of  two  bodies,  as  shown  in  Fig.  3,  are  represented  as: 

a  -  j_b  [ z (b^  - X’t  )  +  Qzi i  )  (1) 

zz  '  12  2 
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Subscripts  1  and  2  refer  to  the  two  bodies  making  the  contact  as  shown  in  Fig.  3. 

R^,  R| ,  R2  and  R2  are  the  principal  radii  of  curvature  at  the  point  of  contact  and 

a  is  the  angle  between  the  corresponding  planes  of  the  principal  curvatures.  "E"  is  the 
modulus  of  elasticity  and  If  is  the  poisson's  ratio. 


Since  the  experiments  were  conducted  for  the  single  chain  of  disks,  the  contact  stress 
field  was  defined  by  only  one  parameter  and  that  is  the  contact  length,  2b.  Fcr  the  two- 
dimensional  plane  stress  problem  the  contact  strains  are  given  by: 


zz 


1  (a  - a  ) 
-  '  zz  *  xx' 


(4) 


Substituting  for  a  and  a from  equation  (1)  and  equation  (2)  and  noting  that  (3  =  0 
(frictionless  case  for  the  single  chain  experiments)  the  contact  strain  is  obtained  as  a 
function  of  the  coordinates  (x  and  z)  and  the  contact  length  (2b) .  Since  the  strain  gages 
are  mounted  normal  to  the  contact  point  the  coordinate  "x"  is  equal  to  zero.  The  contact 
strain  is  plotted  as  a  function  of  the  normal  distance  from  the  conta’ct  for  different  values 
of  (2b)  as  shown  in  Fig.  4.  The  location  of  the  strain  gage  from  the  contact  is  known  (i.e. 
"z"  is  known) .  Also,  the  strain  at  that  location  is  known  from  the  strain  gage  experiment. 
Thus  by  interpolation  the  contact  strain  and  the  contact  length  can  be  obtained. 


Since  the  strain  gages  have  a  finite  size  it  averages  the  strain  over  its  grid  area. 
This  average  strain  is  not  equal  to  the  strain  at  the  grid  geometric  center.  From  the  plot 
of  the  strain  vs  distance  it  can  be  seen  that  steep  gradients  in  the  strain  exist  near  the 
contact  point  and  the  stra-Ln  profiles  peak  around  z  =  1.0  mm.  As  the  distance  normal  to  the 
contact  point  increases  the  strain  gradients  reduce  and  the  strain  profiles  become  fairly 
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constant  after  z  =  6mm.  Thus  it  is  important  that  the  strain  gages  be  placed  beyond  z  =  6 
-r.  so  as  to  minimize  the  averaging  error. 

To  estimate  the  error  induced  due  to  the  averaging  effect  the  percentage  error  is 
calculated  by  considering  a  strain  gage  with  grid  size  (LxL)  positioned  at  a  height  "z" 
above  the  contact  point.  Now,  the  strain  e  at  the  geometric  center  of  the  strain  gage  is 
calculated  from  equation  4.  To  evaluate  the  average  strain  the  strain  gage  grid  is  divided 
into  a  matrix  of  100  x  100  points  and  the  st/ain  at  each  point  is  evaluated.  From  this  the 
average  strain  is  computed.  The  results  are  presented  in  Fig  5,  which  shows  the  error  due 
to  averaging  effect  as  a  function  of  the  gage  grid  length,  L.  The  results  are  plotted  for 
different  heights,  z  above  the  contact  points.  The  error  is  higher  for  strain  gages  closer 
to  the  contact  points  because  of  the  steep  strain  gradients.  Also  the  error  increases  as 
the  strain  gage  size  increases.  For  the  size  and  the  location  of  the  strain  gages  used  in 
this  study  the  error  is  about  2%. 

To  investigate  the  relationship  between  the  wave  velocity  and  the  size  of  the  grains 
a  one  dimensional  assembly  of  disks  as  described  earlier  was  used.  Strain  gages  were  bounded 
on  four  separate  disks  in  the  chain  at  a  height  of  7mm  above  the  contact  point.  The  strain 
gages  were  suitably  connected  to  the  Nicolet  Oscilloscope  through  bridge  amplifiers.  When 
the  explosive  was  fired,  the  oscilloscope  was  triggered.  Strain  as  a  function  of  time  was 
recorded.  The  transit  time  and  the  rise  time  of  the  gage  are  small  (less  than  1  %)  as 
compared  to  the  pulse  duration.  Also  the  frequency  response  of  the  bridge  amplifiers  was 
adequate  for  the  experiments. 

To  study  and  quantify  the  damage  induced  in  the  media  due  to  the  propagation  of  the 
shock  waves  both  the  residual  strength  as  well  as  the  wave  velocity  was  computed  for  the 
single  chain  of  the  disk  assembly  for  repeated  loadings.  The  residual  tensile  strength  was 
measured  by  the  tensile  splitting  tests.  In  these  tests  the  disks  fabricated  from  the  cores 
of  different  grades  of  marble  were  laid  vertically  between  the  loading  plates  of  the  Instron 
machine  in  the  compression  mode.  The  load  was  slowly  increased  at  a  very  slow-rate  until 
the  specimens  failed  by  splitting  across  the  vertical  diameter.  Assuming  linear  condition 
within  the  core  the  approximate  tensile  strength  aT,  was  calculated  using  the  equation: 

2  P 

ct  =  max  ( ->) 

rnDT 

p 

where  max  -  applied  load  at  failure 

D  -  Diameter  of  the  core 
T  -  Thickness  of  the  specimen 

To  estimate  the  residual  strength  of  the  damaged  rock  disks  15  mg  of  Lead  Azide  was 
used  to  generate  the  shock  wave  through  the  assembly  of  disks.  Tensile  splitting  test  was 
performed  on  the  first  two  disks  taken  out  from  the  top  of  the  assembly.  Two  new  disks  were 
placed  on  top  of  the  disk  assembly  to  replace  the  damaged  disks.  The  chain  assembly,  this 
time  was  loaded  twice,  and  again  the  top  two  disks  were  taken  out  and  their  residual  tensile 
strength  was  obtained.  This  procedure  was  continued  until  the  residual  tensile  strength  was 
obtained  for  the  disks  which  had  been  shocked  six  times. 

To  study  the  effect  of  compaction  and  the  damage  of  the  disks  on  the  wave  velocity 
strain  gages  were  bonded  on  the  first,  second,  fourth  and  the  seventh  disks  in  the  chain 
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assembly.  Again,  15mg  of  lead  azide  was  used  to  generate  the  shock  wave  through  the  disk 
assembly.  From  the  peaks  of  the  strain  profile  the  stress  wave  velocity  was  computed  every 
time  the  disk  assembly  was  loaded.  This  procedure  was  repeated  six  times,  keeping  the 
location  of  the  strain  gages  the  same,  and  each  time  average  wave  velocity  and  the  stress 
wave  attenuation  was  computed.  This  average  wave  velocity  versus  the  number  of  times  the 
assembly  was  loaded  was  plotted  to  obtain  the  variation  of  the  wave  velocity  as  a  function 
of  the  load  history. 

RESULTS  AND  DISCUSSIONS 


Wave  Velocity  in  Virgin  Rock  Disk  Assembly 

A  series  of  experiments  were  conducted  with  strain  gages  mounted  7mm  from  the  contact 
point  of  the  disks.  Strain  versus  time  plots  at  the  location  of  the  strain  gages  were 
recorded.  A  typical  strain  gage  output  profile  is  shown  in  Fig.  6.  Position  of  the  wave 
front  versus  time  plots  were  drawn  for  marble  "A",  marble  "C"  and  marble  "D"  disk  assemblies. 
The  average  stress  wave  velocity  was  computed  by  drawing  a  least  sguare  fit  line  passing 
through  the  experimental  points.  The  results  indicate  that  the  average  wave  velocity 
increases  with  the  increasing  grain  size  of  the  marble.  The  average  grain  size 
of  marble  "A"  was  0.32mm  and  the  stress  wave  velocity  was  1210  m;sec.  The  grain  size  of 
marble  "C"  was  higher  and  the  corresponding  stress-wave  velocity  obtained  was  I425m/sec 
while  marble  "D"  had  the  largest  average  grain  size  and  a  stress-wave  velocity  of  1525  m/sec. 
This  trend  is  to  be  expected  as  the  stress  wave  encounters  a  fewer  number  of  grain  boundaries 
or  obstacles  as  it  propagates  through  the  granular  rock  material  having  a  larger  grain  size. 

To  study  the  attenuation  in  the  three  different  grades  of  marble  mentioned  above  the 
normalized  peak  strain  was  plotted  against  the  stress  wave  propagation  distance.  As  seen 
from  Fig.  7  all  the  three  different  grades  of  marble  follow  more  or  less  the  same  pattern 
of  attenuation,  but  the  marble  with  the  smaller  grain  size  shows  a  higher  attenuation.  At 
a  given  contact  the  strain  increased  from  zero  to  a  peak  value  and  then  gradually  decreased 
to  zero.  A  typical  duration  of  the  pulse  at  the  contact  was  75  microseconds.  The  wavelength 
of  the  pulse  was  calculated  by  taking  the  product  of  the  average  wave  speed  and  the  pulse 
duration  at  the  contact.  The  wave  length  obtained  was  95  mm  (around  4  disk  diameter)  .  Due 
to  the  internal  losses  within  the  disk,  energy  spent  in  closing  the  contact  and  some 
frictional  and  reflection  effects,  the  peak  contact  loads  dropped  as  the  distance  from  the 
explosive  loading  'was  increased.  Using  the  values  of  the  normalized  peak  strains  it  is  seen 
that  a  rapid  attenuation  in  peak  strain  occurs  as  the  wave  propagates  through  the  first 
contact.  After  this  the  decay  in  peak  load  is  gradual  with  distance  and  attenuates  to  35i 
of  its  peak  value  in  the  next  nine  disks. 

The  tensile  splitting  tests  were  performed  on  all  the  four  different  grades  of  virgin 
marble  disks.  Table  "A"  lists  the  tensile  failure  strength  of  the  four  disks.  Results 
indicate  that  marble  "C"  had  the  highest  tensile  strength  of  1440  N/mJ  while  marble  "D"  had 
the  lowest  tensile  strength  of  930N/mJ.  Marble  "A"  and  marble  "B"  had  almost  similar  tensile 
strengths  of  1145N/m2  and  1125N/m2  respectively. 


WAVE  PROPAGATION  IN  DAMAGED  ROCK  DISKS 

To  quantify  damage  occurring  in  the  rock  disks,  tensile  splitting  tests  were  carried 
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out  to  obtain  the  residual  tensile  strength.  The  tensile  strength  of  the  first  two  disks 
was  only  determined  since  most  of  the  attenuation  in  stress  waves  occurs  in  these  disks. 
Table  B  lists  the  residual  tensile  strength  of  the  first  two  disks  as  a  function  of  the 
number  of  tines  the  assembly  was  loaded.  It  is  observed  that  after  the  passage  of  the  first 
shock  wave  the  residual  strength  of  the  first  disk  falls  down  to  almost  40%  of  its  value  in 
the  virgin  state.  Additional  drop  of  10  percent  in  the  tensile  strength  occurs  when  the 
second  loading  wave  passes  through  the  disk  assembly.  The  residual  strength  remains  fairly 
constant  as  tne  number  of  stress-wave  loadings  are  increased  fur-ther  and  eventually  the  disk 
fails  by  fracturing  across  the  surface. 

The  tensile  strength  of  the  second  disk  in  the  assembly  drops  down  by  only  10%  after 
the  first  two  stress  waves.  This  happens  because  the  cumulative  damage  occurring  in  the 
second  disk  is  small  as  compared  to  the  first  disk  as  the  peak  stress  value  attenuates  to 
almost  60%  by  the  time  the  stress  wave  propagates  one  disk  diameter.  As  the  number  of  shock 
wave  loading  is  increased  further  the  residual  strength  remains  fairly  constant  until  the 
fourth  loading  wave  after  which  the  tensile  strength  falls  to  72%  of  the  tensile  strength 
value  in  the  virgin  state. 

Next,  the  amount  of  explosive  was  gradually  increased  and  its  effect  on  the  residual 
tensile  strength  studied.  The  results  are  shown  in  Table  C.  When  15  mg  of  lead  Azide  was 
used  the  tensile  strength  of  the  second  disk  falls  down  to  95%  of  its  value  in  the  virgin 
state  after  the  passage  of  one  shock  wave.  When  the  charge  is  increased  to  20mg  the  tensile 
strength  falls  down  to  90%.  As  the  explosive  charge  is  increased  to  50mg  the  tensile  strength 
drops  considerably  to  55%.  In  case  of  20mg,  40mg  and  50mg  of  lead  Azide  the  first  disk 
tailed  by  fracturing  across  the  surface.  It  is  observed  from  all  the  experiments  that  the 
residual  tensile  strength  of  the  disks  approaches  a  limiting  value  of  650N/m2  after  which 
upon  repeated  loading  the  disk  fractures. 

The  average  wave  velocity  in  the  chain  was  calculated  as  a  function  of  repeated 
explosive  loading  to  study  the  effect  of  the  cumulative  damage  on  the  wave  velocity.  Fig. 
3  shows  the  plot  of  the  average  wave  velocity  as  a  function  of  the  number  of  loading  waves. 
The  stress  wave  velocity  is  seen  to  increase  considerably  after  the  passage  of  first  loading 
wave.  The  stress  wave  velocity  increases  to  1400m/sec  from  1080m/sec  which  is  the  velocity 
obtained  by  loading  the  undamaged  disk  chain  assembly.  As  the  number  of  stress  wave  loadings 
is  increased  further  the  wave  velocity  undulates  about  the  terminal  velocity  of  1450m/sec 
and  eventually  falls  down  considerably  as  the  cumulative  damage  increases  and  fractures  the 
disk.  The  increase  in  wave  velocity  with  increased  number  of  stress  wave  loading  is  a  result 
of  the  compacting  process  by  which  the  disks  come  closer  to  each  other  thus  increasing  the 
resultant  sti  fness  of  the  disk  chain  assembly.  This  also  results  in  closing  of  the 
preexisting  voids  and  microcracks  in  the  disks.  The  surface  cracks  impede  the  progress  of 
the  stress-wave  as  it  has  to  travel  around  these  cracks.  Thus  the  compactive  process  re¬ 
duces  the  transit  time  of  propagation  of  stress  waves  resulting  in  an  increase  in  velocity. 


As  the  number  of  shock  wave  loadings  is  increased  damage  occurs  around  the  contact  zone 
due  to  high  contact  stresses  and  large  scale  crumbling  and  pitting  is  observed  near  the 
contact  zone.  Figure  9  shows  the  photomicrographs  of  the  area  near  the  contact  zone  of  the 
disks  after  the  passage  of  the  first  two  stress  waves.  The  stress  wave  velocity  reaches  a 
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terminal  velocity  at  this  stage.  Upon  repeated  loading  the  microcracks  grow  and  engulf  the 
entire  contact  zone.  These  microcracks  combine  together  to  create  open  cracks  on  the  disk 
surface  as  shown  in  photomicrograph  in  Fig.  ie.  This  microcracking  of  the  contact  zone 
lowers  the  stiffness  of  the  assembly  and  the  cracks  impede  the  propagation  of  the  stress  wave 
and  hence  lower  the  stress  wave  velocity. 


CONCLUSIONS 


Dynamic  strain  gage  technique  is  employed  to  study  wave  propagation  and  dynamic  load 
transfer  in  granular  rock  media.  The  granular  rock  media  was  modelled  as  a  one  dimensional 
array  of  circular  disks  fabricated  from  different  grades  of  Vermont  marbles  having  different 
microstructures .  The  dynamic  data  was  analyzed  to  obtain  average  wave  velocity,  stress  wave 
attenuation  for  both  virgin  as  well  as  damaged  rock  media.  The  residual  strength  of  the  rock 
disks  was  estimated  after  the  passage  of  the  stress  wave  using  tensile  splitting  tests.  The 
damage  in  the  disks  was  studied  by  taking  optical  micrographs  of  the  disk  surface.  The 
results  indicate  that: 

1.  A  correlation  exists  between  the  microstructure  of  the  rock  material  and  the  stress 
wave  velocity.  Larger  the  grain  size  of  the  rock  media  higher  is  the  stress-wave 
velocity. 

2.  Rapid  attenuation  in  the  peak  contact  loads,  and  thus  the  amplitude  of  the  stress  wave, 

takes  place  initially  as  the  wave  propagates  through  the  granular  rock  media  and  then 

approaches  a  steady  rate  of  decay.  It  is  observed  that  the  stress  wave  attenuation 
decreases  as  the  average  grain  size  of  the  marble  increases. 

3.  The  stress  wave  velocity  increases  sharply  after  the  passage  of  the  first  shock  wave. 

Upon  repeated  explosive  loading  the  stress-wave  velocity  approaches  a  terminal  velocity 
and  eventually  falls  off  as  the  disk  fails  by  fracturing  across  the  surface. 

4.  The  residual  tensile  strength  of  the  disk  decreases  as  the  number  of  shock  waves  are 

increased.  Further,  a  limiting  residual  strength  value  exists  which  remains  fairly 
constant  upon  repeated  loading  until  the  disk  fractures. 
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Fig  6.  A  typical  strain  profile  obtained  during 
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Fig  9.  Damage  as  a  result  of  stress  wave 
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Fig  10.  Photomicrograph  showing  damage  near  the  contact 
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TABLE  A 


TENSILE  STRENGTH  N/nT)  OBTAINED  FROM  TENSILE  SPLITTING  TEST 


ROCK 

TYPE 

. 

i 

2 

3 

4+ 

average 

STRENGTH 

MARBLE  A 

1130 

1040 

1180 

1180 

in- 

MARBLE  3 

1280 

1020 

1170 

1070 

1125 

MARBLE  C 

U60 

15-0 

1380 

1385 

14-1 

MARBLE  D 

800 

11-0 

1120 

6  70 

93  2 

TABL-  B 


NUMBER  OF  TIMES  THE 

SHOCK  WAVE  PASSED 

THROUGH  THE  DISK 

CHAIN  ASSEMBLY 

(15  mg  OF  LEAD  AZIDE) 

DISK 

NUMBER 

RESIDUAL  TENSILE 
STRENGTH 
(N/m;) 

1 

610 

1 

2 

1160 

1 

610 

2 

2 

980 

1 

640 

3 

2 

it  0 

i 

FRACTURED 

4 

2 

1030 

i 

670 

5 

2 

690 

TABLE  C 


LEAD  AZIDE 

RESIDUAL  TENSILE  STRENGTH 

(rag) 

( N/nr ) 

DISK  =L 

DISK  =2 

15 

610 

1160 

20 

fractured 

1070 

40 

fractured 

780 

50 

fractured 

653 

110 


